R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(-6
+ ,591
+ ,2981.85
+ ,2
+ ,-3
+ ,589
+ ,3080.58
+ ,2
+ ,-2
+ ,584
+ ,3106.22
+ ,2
+ ,-5
+ ,573
+ ,3119.31
+ ,2
+ ,-11
+ ,567
+ ,3061.26
+ ,2
+ ,-11
+ ,569
+ ,3097.31
+ ,2
+ ,-11
+ ,621
+ ,3161.69
+ ,2
+ ,-10
+ ,629
+ ,3257.16
+ ,2
+ ,-14
+ ,628
+ ,3277.01
+ ,2
+ ,-8
+ ,612
+ ,3295.32
+ ,2
+ ,-9
+ ,595
+ ,3363.99
+ ,2
+ ,-5
+ ,597
+ ,3494.17
+ ,2.21
+ ,-1
+ ,593
+ ,3667.03
+ ,2.25
+ ,-2
+ ,590
+ ,3813.06
+ ,2.25
+ ,-5
+ ,580
+ ,3917.96
+ ,2.45
+ ,-4
+ ,574
+ ,3895.51
+ ,2.5
+ ,-6
+ ,573
+ ,3801.06
+ ,2.5
+ ,-2
+ ,573
+ ,3570.12
+ ,2.64
+ ,-2
+ ,620
+ ,3701.61
+ ,2.75
+ ,-2
+ ,626
+ ,3862.27
+ ,2.93
+ ,-2
+ ,620
+ ,3970.1
+ ,3
+ ,2
+ ,588
+ ,4138.52
+ ,3.17
+ ,1
+ ,566
+ ,4199.75
+ ,3.25
+ ,-8
+ ,557
+ ,4290.89
+ ,3.39
+ ,-1
+ ,561
+ ,4443.91
+ ,3.5
+ ,1
+ ,549
+ ,4502.64
+ ,3.5
+ ,-1
+ ,532
+ ,4356.98
+ ,3.65
+ ,2
+ ,526
+ ,4591.27
+ ,3.75
+ ,2
+ ,511
+ ,4696.96
+ ,3.75
+ ,1
+ ,499
+ ,4621.4
+ ,3.9
+ ,-1
+ ,555
+ ,4562.84
+ ,4
+ ,-2
+ ,565
+ ,4202.52
+ ,4
+ ,-2
+ ,542
+ ,4296.49
+ ,4
+ ,-1
+ ,527
+ ,4435.23
+ ,4
+ ,-8
+ ,510
+ ,4105.18
+ ,4
+ ,-4
+ ,514
+ ,4116.68
+ ,4
+ ,-6
+ ,517
+ ,3844.49
+ ,4
+ ,-3
+ ,508
+ ,3720.98
+ ,4
+ ,-3
+ ,493
+ ,3674.4
+ ,4
+ ,-7
+ ,490
+ ,3857.62
+ ,4
+ ,-9
+ ,469
+ ,3801.06
+ ,4
+ ,-11
+ ,478
+ ,3504.37
+ ,4
+ ,-13
+ ,528
+ ,3032.6
+ ,4.18
+ ,-11
+ ,534
+ ,3047.03
+ ,4.25
+ ,-9
+ ,518
+ ,2962.34
+ ,4.25
+ ,-17
+ ,506
+ ,2197.82
+ ,3.97
+ ,-22
+ ,502
+ ,2014.45
+ ,3.42
+ ,-25
+ ,516
+ ,1862.83
+ ,2.75
+ ,-20
+ ,528
+ ,1905.41
+ ,2.31
+ ,-24
+ ,533
+ ,1810.99
+ ,2
+ ,-24
+ ,536
+ ,1670.07
+ ,1.66
+ ,-22
+ ,537
+ ,1864.44
+ ,1.31
+ ,-19
+ ,524
+ ,2052.02
+ ,1.09
+ ,-18
+ ,536
+ ,2029.6
+ ,1
+ ,-17
+ ,587
+ ,2070.83
+ ,1
+ ,-11
+ ,597
+ ,2293.41
+ ,1
+ ,-11
+ ,581
+ ,2443.27
+ ,1
+ ,-12
+ ,564
+ ,2513.17
+ ,1
+ ,-10
+ ,558
+ ,2466.92
+ ,1
+ ,-15
+ ,575
+ ,2502.66
+ ,1
+ ,-15
+ ,580
+ ,2539.91
+ ,1
+ ,-15
+ ,575
+ ,2482.6
+ ,1
+ ,-13
+ ,563
+ ,2626.15
+ ,1
+ ,-8
+ ,552
+ ,2656.32
+ ,1
+ ,-13
+ ,537
+ ,2446.66
+ ,1
+ ,-9
+ ,545
+ ,2467.38
+ ,1
+ ,-7
+ ,601
+ ,2462.32
+ ,1
+ ,-4
+ ,604
+ ,2504.58
+ ,1
+ ,-4
+ ,586
+ ,2579.39
+ ,1
+ ,-2
+ ,564
+ ,2649.24
+ ,1
+ ,0
+ ,549
+ ,2636.87
+ ,1
+ ,-2
+ ,551
+ ,2613.94
+ ,1
+ ,-3
+ ,556
+ ,2634.01
+ ,1
+ ,1
+ ,548
+ ,2711.94
+ ,1
+ ,-2
+ ,540
+ ,2646.43
+ ,1
+ ,-1
+ ,531
+ ,2717.79
+ ,1.14
+ ,1
+ ,521
+ ,2701.54
+ ,1.25
+ ,-3
+ ,519
+ ,2572.98
+ ,1.25
+ ,-4
+ ,572
+ ,2488.92
+ ,1.4
+ ,-9
+ ,581
+ ,2204.91
+ ,1.5
+ ,-9
+ ,563
+ ,2123.99
+ ,1.5
+ ,-7
+ ,548
+ ,2149.1
+ ,1.5)
+ ,dim=c(4
+ ,82)
+ ,dimnames=list(c('Consumentenvertrouwen'
+ ,'Werkloosheid'
+ ,'BEL20'
+ ,'Rentevoet')
+ ,1:82))
> y <- array(NA,dim=c(4,82),dimnames=list(c('Consumentenvertrouwen','Werkloosheid','BEL20','Rentevoet'),1:82))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> library(lattice)
> library(lmtest)
Loading required package: zoo
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Consumentenvertrouwen Werkloosheid BEL20 Rentevoet
1 -6 591 2981.85 2.00
2 -3 589 3080.58 2.00
3 -2 584 3106.22 2.00
4 -5 573 3119.31 2.00
5 -11 567 3061.26 2.00
6 -11 569 3097.31 2.00
7 -11 621 3161.69 2.00
8 -10 629 3257.16 2.00
9 -14 628 3277.01 2.00
10 -8 612 3295.32 2.00
11 -9 595 3363.99 2.00
12 -5 597 3494.17 2.21
13 -1 593 3667.03 2.25
14 -2 590 3813.06 2.25
15 -5 580 3917.96 2.45
16 -4 574 3895.51 2.50
17 -6 573 3801.06 2.50
18 -2 573 3570.12 2.64
19 -2 620 3701.61 2.75
20 -2 626 3862.27 2.93
21 -2 620 3970.10 3.00
22 2 588 4138.52 3.17
23 1 566 4199.75 3.25
24 -8 557 4290.89 3.39
25 -1 561 4443.91 3.50
26 1 549 4502.64 3.50
27 -1 532 4356.98 3.65
28 2 526 4591.27 3.75
29 2 511 4696.96 3.75
30 1 499 4621.40 3.90
31 -1 555 4562.84 4.00
32 -2 565 4202.52 4.00
33 -2 542 4296.49 4.00
34 -1 527 4435.23 4.00
35 -8 510 4105.18 4.00
36 -4 514 4116.68 4.00
37 -6 517 3844.49 4.00
38 -3 508 3720.98 4.00
39 -3 493 3674.40 4.00
40 -7 490 3857.62 4.00
41 -9 469 3801.06 4.00
42 -11 478 3504.37 4.00
43 -13 528 3032.60 4.18
44 -11 534 3047.03 4.25
45 -9 518 2962.34 4.25
46 -17 506 2197.82 3.97
47 -22 502 2014.45 3.42
48 -25 516 1862.83 2.75
49 -20 528 1905.41 2.31
50 -24 533 1810.99 2.00
51 -24 536 1670.07 1.66
52 -22 537 1864.44 1.31
53 -19 524 2052.02 1.09
54 -18 536 2029.60 1.00
55 -17 587 2070.83 1.00
56 -11 597 2293.41 1.00
57 -11 581 2443.27 1.00
58 -12 564 2513.17 1.00
59 -10 558 2466.92 1.00
60 -15 575 2502.66 1.00
61 -15 580 2539.91 1.00
62 -15 575 2482.60 1.00
63 -13 563 2626.15 1.00
64 -8 552 2656.32 1.00
65 -13 537 2446.66 1.00
66 -9 545 2467.38 1.00
67 -7 601 2462.32 1.00
68 -4 604 2504.58 1.00
69 -4 586 2579.39 1.00
70 -2 564 2649.24 1.00
71 0 549 2636.87 1.00
72 -2 551 2613.94 1.00
73 -3 556 2634.01 1.00
74 1 548 2711.94 1.00
75 -2 540 2646.43 1.00
76 -1 531 2717.79 1.14
77 1 521 2701.54 1.25
78 -3 519 2572.98 1.25
79 -4 572 2488.92 1.40
80 -9 581 2204.91 1.50
81 -9 563 2123.99 1.50
82 -7 548 2149.10 1.50
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Werkloosheid BEL20 Rentevoet
-15.206252 -0.023441 0.009417 -3.731899
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-7.4913 -3.6352 -0.4637 3.9753 7.6448
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.521e+01 8.764e+00 -1.735 0.0867 .
Werkloosheid -2.344e-02 1.591e-02 -1.474 0.1446
BEL20 9.417e-03 9.077e-04 10.374 2.42e-16 ***
Rentevoet -3.732e+00 7.062e-01 -5.284 1.11e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.241 on 78 degrees of freedom
Multiple R-squared: 0.6257, Adjusted R-squared: 0.6113
F-statistic: 43.47 on 3 and 78 DF, p-value: < 2.2e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.817936384 0.3641272311 0.1820636156
[2,] 0.701536271 0.5969274574 0.2984637287
[3,] 0.653523841 0.6929523175 0.3464761587
[4,] 0.588182316 0.8236353688 0.4118176844
[5,] 0.488394808 0.9767896166 0.5116051917
[6,] 0.375913352 0.7518267033 0.6240866483
[7,] 0.312054787 0.6241095730 0.6879452135
[8,] 0.229594523 0.4591890455 0.7704054772
[9,] 0.307123447 0.6142468934 0.6928765533
[10,] 0.260478123 0.5209562454 0.7395218773
[11,] 0.231585174 0.4631703477 0.7684148262
[12,] 0.180785367 0.3615707348 0.8192146326
[13,] 0.138274973 0.2765499464 0.8617250268
[14,] 0.095916946 0.1918338917 0.9040830541
[15,] 0.065197341 0.1303946824 0.9348026588
[16,] 0.044153265 0.0883065294 0.9558467353
[17,] 0.031083367 0.0621667339 0.9689166330
[18,] 0.192775046 0.3855500924 0.8072249538
[19,] 0.154098277 0.3081965545 0.8459017227
[20,] 0.114933480 0.2298669605 0.8850665198
[21,] 0.088063616 0.1761272327 0.9119363837
[22,] 0.062619594 0.1252391887 0.9373804056
[23,] 0.044089698 0.0881793963 0.9559103019
[24,] 0.031083303 0.0621666055 0.9689166973
[25,] 0.025041810 0.0500836196 0.9749581902
[26,] 0.017249085 0.0344981707 0.9827509147
[27,] 0.012108197 0.0242163948 0.9878918026
[28,] 0.008438291 0.0168765824 0.9915617088
[29,] 0.016491884 0.0329837681 0.9835081160
[30,] 0.012523840 0.0250476796 0.9874761602
[31,] 0.009078374 0.0181567485 0.9909216257
[32,] 0.006924962 0.0138499245 0.9930750377
[33,] 0.004996571 0.0099931414 0.9950034293
[34,] 0.004740162 0.0094803240 0.9952598380
[35,] 0.008314818 0.0166296365 0.9916851818
[36,] 0.018731706 0.0374634117 0.9812682942
[37,] 0.017848847 0.0356976934 0.9821511533
[38,] 0.017162644 0.0343252879 0.9828373560
[39,] 0.020370835 0.0407416691 0.9796291654
[40,] 0.013807327 0.0276146531 0.9861926734
[41,] 0.016351254 0.0327025073 0.9836487463
[42,] 0.039679739 0.0793594778 0.9603202611
[43,] 0.045475525 0.0909510492 0.9545244754
[44,] 0.119157269 0.2383145388 0.8808427306
[45,] 0.118123417 0.2362468333 0.8818765834
[46,] 0.101139317 0.2022786338 0.8988606831
[47,] 0.077206736 0.1544134720 0.9227932640
[48,] 0.058012889 0.1160257784 0.9419871108
[49,] 0.047307832 0.0946156649 0.9526921676
[50,] 0.055308087 0.1106161731 0.9446919135
[51,] 0.042820912 0.0856418243 0.9571790878
[52,] 0.035902674 0.0718053480 0.9640973260
[53,] 0.029162306 0.0583246125 0.9708376937
[54,] 0.043820822 0.0876416439 0.9561791780
[55,] 0.115748788 0.2314975754 0.8842512123
[56,] 0.235707704 0.4714154079 0.7642922960
[57,] 0.756142707 0.4877145863 0.2438572932
[58,] 0.915021095 0.1699578104 0.0849789052
[59,] 0.991140001 0.0177199976 0.0088599988
[60,] 0.999466076 0.0010678488 0.0005339244
[61,] 0.999505406 0.0009891885 0.0004945942
[62,] 0.999300735 0.0013985292 0.0006992646
[63,] 0.998657871 0.0026842571 0.0013421285
[64,] 0.997366889 0.0052662212 0.0026331106
[65,] 0.997815404 0.0043691928 0.0021845964
[66,] 0.994358613 0.0112827734 0.0056413867
[67,] 0.987483578 0.0250328433 0.0125164217
[68,] 0.990602747 0.0187945068 0.0093972534
[69,] 0.983485999 0.0330280022 0.0165140011
> postscript(file="/var/www/rcomp/tmp/1xi271321886398.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/2zyb61321886398.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/3nnq31321886398.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/4aube1321886398.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/5f6at1321886398.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 82
Frequency = 1
1 2 3 4 5 6
2.44506264 4.46848098 5.10983308 1.72871512 -3.86530056 -4.15788555
7 8 9 10 11 12
-3.54517438 -3.25664441 -7.46700488 -2.01448418 -4.05962381 -0.45489276
13 14 15 16 17 18
1.97286774 -0.47255951 -3.94839260 -2.69104359 -3.82508897 2.87204183
19 20 21 22 23 24
3.14610839 2.44563055 1.55082549 3.84918441 2.05544924 -7.49128426
25 26 27 28 29 30
-1.42793505 -0.26226709 -0.72936632 0.29696487 -1.04989383 -1.05988854
31 32 33 34 35 36
-0.82254772 1.80484796 0.38082071 -0.27725596 -4.56781681 -0.58234205
37 38 39 40 41 42
0.05108093 4.00315051 4.09015444 -1.70547542 -3.66514217 -2.66036483
43 44 45 46 47 48
1.62590476 3.89190456 6.31433300 4.18726883 -1.23232259 -4.97677383
49 50 51 52 53 54
-1.73847093 -5.88903924 -5.76057626 -6.87360009 -6.76571752 -5.60917238
55 56 57 58 59 60
-3.80190874 0.33656246 -1.44966794 -3.50638996 -1.21152128 -6.14956677
61 62 63 64 65 66
-6.38312759 -5.96067018 -5.59371653 -1.13566968 -4.51300988 -0.52059056
67 68 69 70 71 72
2.83977336 5.51215276 4.38575398 5.21229600 6.97715873 5.23996362
73 74 75 76 77 78
4.16817965 7.24681471 4.67616400 5.31569045 7.64480515 4.80851796
79 80 81 82
6.40225243 4.66081740 5.00086262 6.41279193
> postscript(file="/var/www/rcomp/tmp/6wh0u1321886398.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 82
Frequency = 1
lag(myerror, k = 1) myerror
0 2.44506264 NA
1 4.46848098 2.44506264
2 5.10983308 4.46848098
3 1.72871512 5.10983308
4 -3.86530056 1.72871512
5 -4.15788555 -3.86530056
6 -3.54517438 -4.15788555
7 -3.25664441 -3.54517438
8 -7.46700488 -3.25664441
9 -2.01448418 -7.46700488
10 -4.05962381 -2.01448418
11 -0.45489276 -4.05962381
12 1.97286774 -0.45489276
13 -0.47255951 1.97286774
14 -3.94839260 -0.47255951
15 -2.69104359 -3.94839260
16 -3.82508897 -2.69104359
17 2.87204183 -3.82508897
18 3.14610839 2.87204183
19 2.44563055 3.14610839
20 1.55082549 2.44563055
21 3.84918441 1.55082549
22 2.05544924 3.84918441
23 -7.49128426 2.05544924
24 -1.42793505 -7.49128426
25 -0.26226709 -1.42793505
26 -0.72936632 -0.26226709
27 0.29696487 -0.72936632
28 -1.04989383 0.29696487
29 -1.05988854 -1.04989383
30 -0.82254772 -1.05988854
31 1.80484796 -0.82254772
32 0.38082071 1.80484796
33 -0.27725596 0.38082071
34 -4.56781681 -0.27725596
35 -0.58234205 -4.56781681
36 0.05108093 -0.58234205
37 4.00315051 0.05108093
38 4.09015444 4.00315051
39 -1.70547542 4.09015444
40 -3.66514217 -1.70547542
41 -2.66036483 -3.66514217
42 1.62590476 -2.66036483
43 3.89190456 1.62590476
44 6.31433300 3.89190456
45 4.18726883 6.31433300
46 -1.23232259 4.18726883
47 -4.97677383 -1.23232259
48 -1.73847093 -4.97677383
49 -5.88903924 -1.73847093
50 -5.76057626 -5.88903924
51 -6.87360009 -5.76057626
52 -6.76571752 -6.87360009
53 -5.60917238 -6.76571752
54 -3.80190874 -5.60917238
55 0.33656246 -3.80190874
56 -1.44966794 0.33656246
57 -3.50638996 -1.44966794
58 -1.21152128 -3.50638996
59 -6.14956677 -1.21152128
60 -6.38312759 -6.14956677
61 -5.96067018 -6.38312759
62 -5.59371653 -5.96067018
63 -1.13566968 -5.59371653
64 -4.51300988 -1.13566968
65 -0.52059056 -4.51300988
66 2.83977336 -0.52059056
67 5.51215276 2.83977336
68 4.38575398 5.51215276
69 5.21229600 4.38575398
70 6.97715873 5.21229600
71 5.23996362 6.97715873
72 4.16817965 5.23996362
73 7.24681471 4.16817965
74 4.67616400 7.24681471
75 5.31569045 4.67616400
76 7.64480515 5.31569045
77 4.80851796 7.64480515
78 6.40225243 4.80851796
79 4.66081740 6.40225243
80 5.00086262 4.66081740
81 6.41279193 5.00086262
82 NA 6.41279193
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 4.46848098 2.44506264
[2,] 5.10983308 4.46848098
[3,] 1.72871512 5.10983308
[4,] -3.86530056 1.72871512
[5,] -4.15788555 -3.86530056
[6,] -3.54517438 -4.15788555
[7,] -3.25664441 -3.54517438
[8,] -7.46700488 -3.25664441
[9,] -2.01448418 -7.46700488
[10,] -4.05962381 -2.01448418
[11,] -0.45489276 -4.05962381
[12,] 1.97286774 -0.45489276
[13,] -0.47255951 1.97286774
[14,] -3.94839260 -0.47255951
[15,] -2.69104359 -3.94839260
[16,] -3.82508897 -2.69104359
[17,] 2.87204183 -3.82508897
[18,] 3.14610839 2.87204183
[19,] 2.44563055 3.14610839
[20,] 1.55082549 2.44563055
[21,] 3.84918441 1.55082549
[22,] 2.05544924 3.84918441
[23,] -7.49128426 2.05544924
[24,] -1.42793505 -7.49128426
[25,] -0.26226709 -1.42793505
[26,] -0.72936632 -0.26226709
[27,] 0.29696487 -0.72936632
[28,] -1.04989383 0.29696487
[29,] -1.05988854 -1.04989383
[30,] -0.82254772 -1.05988854
[31,] 1.80484796 -0.82254772
[32,] 0.38082071 1.80484796
[33,] -0.27725596 0.38082071
[34,] -4.56781681 -0.27725596
[35,] -0.58234205 -4.56781681
[36,] 0.05108093 -0.58234205
[37,] 4.00315051 0.05108093
[38,] 4.09015444 4.00315051
[39,] -1.70547542 4.09015444
[40,] -3.66514217 -1.70547542
[41,] -2.66036483 -3.66514217
[42,] 1.62590476 -2.66036483
[43,] 3.89190456 1.62590476
[44,] 6.31433300 3.89190456
[45,] 4.18726883 6.31433300
[46,] -1.23232259 4.18726883
[47,] -4.97677383 -1.23232259
[48,] -1.73847093 -4.97677383
[49,] -5.88903924 -1.73847093
[50,] -5.76057626 -5.88903924
[51,] -6.87360009 -5.76057626
[52,] -6.76571752 -6.87360009
[53,] -5.60917238 -6.76571752
[54,] -3.80190874 -5.60917238
[55,] 0.33656246 -3.80190874
[56,] -1.44966794 0.33656246
[57,] -3.50638996 -1.44966794
[58,] -1.21152128 -3.50638996
[59,] -6.14956677 -1.21152128
[60,] -6.38312759 -6.14956677
[61,] -5.96067018 -6.38312759
[62,] -5.59371653 -5.96067018
[63,] -1.13566968 -5.59371653
[64,] -4.51300988 -1.13566968
[65,] -0.52059056 -4.51300988
[66,] 2.83977336 -0.52059056
[67,] 5.51215276 2.83977336
[68,] 4.38575398 5.51215276
[69,] 5.21229600 4.38575398
[70,] 6.97715873 5.21229600
[71,] 5.23996362 6.97715873
[72,] 4.16817965 5.23996362
[73,] 7.24681471 4.16817965
[74,] 4.67616400 7.24681471
[75,] 5.31569045 4.67616400
[76,] 7.64480515 5.31569045
[77,] 4.80851796 7.64480515
[78,] 6.40225243 4.80851796
[79,] 4.66081740 6.40225243
[80,] 5.00086262 4.66081740
[81,] 6.41279193 5.00086262
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 4.46848098 2.44506264
2 5.10983308 4.46848098
3 1.72871512 5.10983308
4 -3.86530056 1.72871512
5 -4.15788555 -3.86530056
6 -3.54517438 -4.15788555
7 -3.25664441 -3.54517438
8 -7.46700488 -3.25664441
9 -2.01448418 -7.46700488
10 -4.05962381 -2.01448418
11 -0.45489276 -4.05962381
12 1.97286774 -0.45489276
13 -0.47255951 1.97286774
14 -3.94839260 -0.47255951
15 -2.69104359 -3.94839260
16 -3.82508897 -2.69104359
17 2.87204183 -3.82508897
18 3.14610839 2.87204183
19 2.44563055 3.14610839
20 1.55082549 2.44563055
21 3.84918441 1.55082549
22 2.05544924 3.84918441
23 -7.49128426 2.05544924
24 -1.42793505 -7.49128426
25 -0.26226709 -1.42793505
26 -0.72936632 -0.26226709
27 0.29696487 -0.72936632
28 -1.04989383 0.29696487
29 -1.05988854 -1.04989383
30 -0.82254772 -1.05988854
31 1.80484796 -0.82254772
32 0.38082071 1.80484796
33 -0.27725596 0.38082071
34 -4.56781681 -0.27725596
35 -0.58234205 -4.56781681
36 0.05108093 -0.58234205
37 4.00315051 0.05108093
38 4.09015444 4.00315051
39 -1.70547542 4.09015444
40 -3.66514217 -1.70547542
41 -2.66036483 -3.66514217
42 1.62590476 -2.66036483
43 3.89190456 1.62590476
44 6.31433300 3.89190456
45 4.18726883 6.31433300
46 -1.23232259 4.18726883
47 -4.97677383 -1.23232259
48 -1.73847093 -4.97677383
49 -5.88903924 -1.73847093
50 -5.76057626 -5.88903924
51 -6.87360009 -5.76057626
52 -6.76571752 -6.87360009
53 -5.60917238 -6.76571752
54 -3.80190874 -5.60917238
55 0.33656246 -3.80190874
56 -1.44966794 0.33656246
57 -3.50638996 -1.44966794
58 -1.21152128 -3.50638996
59 -6.14956677 -1.21152128
60 -6.38312759 -6.14956677
61 -5.96067018 -6.38312759
62 -5.59371653 -5.96067018
63 -1.13566968 -5.59371653
64 -4.51300988 -1.13566968
65 -0.52059056 -4.51300988
66 2.83977336 -0.52059056
67 5.51215276 2.83977336
68 4.38575398 5.51215276
69 5.21229600 4.38575398
70 6.97715873 5.21229600
71 5.23996362 6.97715873
72 4.16817965 5.23996362
73 7.24681471 4.16817965
74 4.67616400 7.24681471
75 5.31569045 4.67616400
76 7.64480515 5.31569045
77 4.80851796 7.64480515
78 6.40225243 4.80851796
79 4.66081740 6.40225243
80 5.00086262 4.66081740
81 6.41279193 5.00086262
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/77or31321886398.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/8estc1321886398.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/9bvjb1321886398.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/rcomp/tmp/10uvhi1321886398.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/11i7lt1321886398.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/122x8h1321886398.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/134ia01321886398.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/140av71321886398.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/rcomp/tmp/15zmbd1321886398.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/rcomp/tmp/16f4qq1321886398.tab")
+ }
>
> try(system("convert tmp/1xi271321886398.ps tmp/1xi271321886398.png",intern=TRUE))
character(0)
> try(system("convert tmp/2zyb61321886398.ps tmp/2zyb61321886398.png",intern=TRUE))
character(0)
> try(system("convert tmp/3nnq31321886398.ps tmp/3nnq31321886398.png",intern=TRUE))
character(0)
> try(system("convert tmp/4aube1321886398.ps tmp/4aube1321886398.png",intern=TRUE))
character(0)
> try(system("convert tmp/5f6at1321886398.ps tmp/5f6at1321886398.png",intern=TRUE))
character(0)
> try(system("convert tmp/6wh0u1321886398.ps tmp/6wh0u1321886398.png",intern=TRUE))
character(0)
> try(system("convert tmp/77or31321886398.ps tmp/77or31321886398.png",intern=TRUE))
character(0)
> try(system("convert tmp/8estc1321886398.ps tmp/8estc1321886398.png",intern=TRUE))
character(0)
> try(system("convert tmp/9bvjb1321886398.ps tmp/9bvjb1321886398.png",intern=TRUE))
character(0)
> try(system("convert tmp/10uvhi1321886398.ps tmp/10uvhi1321886398.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.320 0.300 4.704