R version 2.13.0 (2011-04-13)
Copyright (C) 2011 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(63.46
+ ,60635408600
+ ,7175.79
+ ,8450000
+ ,64.98
+ ,64463389370
+ ,7533.41
+ ,8557000
+ ,65.48
+ ,66750859382
+ ,7749.11
+ ,8614000
+ ,66.35
+ ,69032584255.95
+ ,7990.47
+ ,8639369
+ ,66.80
+ ,72841303906.58
+ ,8393.42
+ ,8678386
+ ,68.00
+ ,72838685607.53
+ ,8343.11
+ ,8730405
+ ,68.37
+ ,75887977453.76
+ ,8645.37
+ ,8777873
+ ,68.63
+ ,78936158549.66
+ ,8950.31
+ ,8819380
+ ,68.58
+ ,81986612544.38
+ ,9244.73
+ ,8868475
+ ,68.87
+ ,85038924238.38
+ ,9529.4
+ ,8923845
+ ,69.24
+ ,87328862789.67
+ ,9714.96
+ ,8989111
+ ,69.93
+ ,85794984680.72
+ ,9477.27
+ ,9052707
+ ,70.33
+ ,88082838620.17
+ ,9675.47
+ ,9103729
+ ,69.65
+ ,91885510657.4
+ ,10076.6
+ ,9118700
+ ,70.52
+ ,96355573323.8
+ ,10512.51
+ ,9165800
+ ,70.25
+ ,101321342608.8
+ ,10991.21
+ ,9218400
+ ,70.06
+ ,105791395836.8
+ ,11396.13
+ ,9283100
+ ,70.73
+ ,113241494103
+ ,12089.41
+ ,9367000
+ ,70.58
+ ,117215849652.8
+ ,12406.29
+ ,9448100
+ ,70.65
+ ,120940898880.6
+ ,12720.18
+ ,9507800
+ ,70.94
+ ,125658810316.5
+ ,13149.04
+ ,9556500
+ ,70.63
+ ,130873880200.8
+ ,13647.2
+ ,9589800
+ ,70.71
+ ,139583536623.4
+ ,14520.74
+ ,9612700
+ ,70.97
+ ,148226520849
+ ,15379.71
+ ,9637800
+ ,71.10
+ ,153789461350
+ ,15899.66
+ ,9672500
+ ,71.44
+ ,161871513310.4
+ ,16672.14
+ ,9709100
+ ,71.65
+ ,171781295610.4
+ ,17639.58
+ ,9738400
+ ,72.01
+ ,178996585758.2
+ ,18325.17
+ ,9767800
+ ,72.00
+ ,176620960002
+ ,18032.12
+ ,9794800
+ ,72.15
+ ,186604680395
+ ,19019.94
+ ,9811000
+ ,72.80
+ ,187772917033.2
+ ,19117.97
+ ,9821800
+ ,72.75
+ ,193109744878.6
+ ,19645.54
+ ,9829700
+ ,73.24
+ ,197630528464
+ ,20090.12
+ ,9837200
+ ,73.29
+ ,206483683756.4
+ ,20969.62
+ ,9846800
+ ,73.70
+ ,203906590621.6
+ ,20696.13
+ ,9852400
+ ,73.93
+ ,206783719069.34
+ ,20979.85
+ ,9856303
+ ,73.93
+ ,206759082201.76
+ ,20979.01
+ ,9855520
+ ,74.43
+ ,211878483671.4
+ ,21498.94
+ ,9855300
+ ,74.54
+ ,214009149959
+ ,21708.75
+ ,9858200
+ ,74.74
+ ,217203709135.4
+ ,22024.75
+ ,9861800
+ ,75.35
+ ,222331811922.6
+ ,22525.56
+ ,9870200
+ ,75.66
+ ,232825724880
+ ,23555.82
+ ,9884000
+ ,75.73
+ ,241180274038.7
+ ,24269.23
+ ,9937697
+ ,76.14
+ ,248494123822.1
+ ,24925.91
+ ,9969310
+ ,76.30
+ ,253049202253.08
+ ,25293.57
+ ,10004487
+ ,76.46
+ ,256922518700.16
+ ,25575.57
+ ,10045622
+ ,76.47
+ ,254451243638.75
+ ,25229.6
+ ,10085426
+ ,76.82
+ ,262662316800.94
+ ,25947.3
+ ,10122914
+ ,76.97
+ ,268926171539.67
+ ,26480.95
+ ,10155459
+ ,77.31
+ ,272037080259.13
+ ,26725.5
+ ,10178934
+ ,77.53
+ ,281118338865.04
+ ,27561.2
+ ,10199787
+ ,77.62
+ ,286389613939.39
+ ,28030.61
+ ,10217030
+ ,77.80
+ ,296190694318.91
+ ,28937.15
+ ,10235655
+ ,77.91
+ ,307294826961.69
+ ,29940.2
+ ,10263618
+ ,78.22
+ ,309697986635.6
+ ,30092.08
+ ,10291679
+ ,78.32
+ ,314369521231.48
+ ,30485.88
+ ,10311970
+ ,78.47
+ ,317533640477.42
+ ,30736.53
+ ,10330824
+ ,79.12
+ ,327005697520.69
+ ,31600.02
+ ,10348276
+ ,79.21
+ ,332458473876
+ ,32077
+ ,10364388
+ ,79.55
+ ,339794119726.27
+ ,32738.41
+ ,10379067)
+ ,dim=c(4
+ ,60)
+ ,dimnames=list(c('Levensverwachting'
+ ,'GDP'
+ ,'Inkomen'
+ ,'Populatie')
+ ,1:60))
> y <- array(NA,dim=c(4,60),dimnames=list(c('Levensverwachting','GDP','Inkomen','Populatie'),1:60))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> 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
Levensverwachting GDP Inkomen Populatie t
1 63.46 60635408600 7175.79 8450000 1
2 64.98 64463389370 7533.41 8557000 2
3 65.48 66750859382 7749.11 8614000 3
4 66.35 69032584256 7990.47 8639369 4
5 66.80 72841303907 8393.42 8678386 5
6 68.00 72838685608 8343.11 8730405 6
7 68.37 75887977454 8645.37 8777873 7
8 68.63 78936158550 8950.31 8819380 8
9 68.58 81986612544 9244.73 8868475 9
10 68.87 85038924238 9529.40 8923845 10
11 69.24 87328862790 9714.96 8989111 11
12 69.93 85794984681 9477.27 9052707 12
13 70.33 88082838620 9675.47 9103729 13
14 69.65 91885510657 10076.60 9118700 14
15 70.52 96355573324 10512.51 9165800 15
16 70.25 101321342609 10991.21 9218400 16
17 70.06 105791395837 11396.13 9283100 17
18 70.73 113241494103 12089.41 9367000 18
19 70.58 117215849653 12406.29 9448100 19
20 70.65 120940898881 12720.18 9507800 20
21 70.94 125658810316 13149.04 9556500 21
22 70.63 130873880201 13647.20 9589800 22
23 70.71 139583536623 14520.74 9612700 23
24 70.97 148226520849 15379.71 9637800 24
25 71.10 153789461350 15899.66 9672500 25
26 71.44 161871513310 16672.14 9709100 26
27 71.65 171781295610 17639.58 9738400 27
28 72.01 178996585758 18325.17 9767800 28
29 72.00 176620960002 18032.12 9794800 29
30 72.15 186604680395 19019.94 9811000 30
31 72.80 187772917033 19117.97 9821800 31
32 72.75 193109744879 19645.54 9829700 32
33 73.24 197630528464 20090.12 9837200 33
34 73.29 206483683756 20969.62 9846800 34
35 73.70 203906590622 20696.13 9852400 35
36 73.93 206783719069 20979.85 9856303 36
37 73.93 206759082202 20979.01 9855520 37
38 74.43 211878483671 21498.94 9855300 38
39 74.54 214009149959 21708.75 9858200 39
40 74.74 217203709135 22024.75 9861800 40
41 75.35 222331811923 22525.56 9870200 41
42 75.66 232825724880 23555.82 9884000 42
43 75.73 241180274039 24269.23 9937697 43
44 76.14 248494123822 24925.91 9969310 44
45 76.30 253049202253 25293.57 10004487 45
46 76.46 256922518700 25575.57 10045622 46
47 76.47 254451243639 25229.60 10085426 47
48 76.82 262662316801 25947.30 10122914 48
49 76.97 268926171540 26480.95 10155459 49
50 77.31 272037080259 26725.50 10178934 50
51 77.53 281118338865 27561.20 10199787 51
52 77.62 286389613939 28030.61 10217030 52
53 77.80 296190694319 28937.15 10235655 53
54 77.91 307294826962 29940.20 10263618 54
55 78.22 309697986636 30092.08 10291679 55
56 78.32 314369521231 30485.88 10311970 56
57 78.47 317533640477 30736.53 10330824 57
58 79.12 327005697521 31600.02 10348276 58
59 79.21 332458473876 32077.00 10364388 59
60 79.55 339794119726 32738.41 10379067 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) GDP Inkomen Populatie t
7.922e+01 2.396e-11 -8.096e-04 -1.151e-06 5.036e-01
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.17497 -0.24768 0.01103 0.30998 0.91530
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.922e+01 5.197e+00 15.244 < 2e-16 ***
GDP 2.396e-11 3.874e-11 0.618 0.5389
Inkomen -8.096e-04 4.326e-04 -1.871 0.0666 .
Populatie -1.151e-06 5.762e-07 -1.998 0.0507 .
t 5.036e-01 5.403e-02 9.320 6.54e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5237 on 55 degrees of freedom
Multiple R-squared: 0.9836, Adjusted R-squared: 0.9825
F-statistic: 827.2 on 4 and 55 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.4016646 8.033292e-01 5.983354e-01
[2,] 0.3443147 6.886295e-01 6.556853e-01
[3,] 0.3227120 6.454240e-01 6.772880e-01
[4,] 0.3065381 6.130762e-01 6.934619e-01
[5,] 0.4731168 9.462336e-01 5.268832e-01
[6,] 0.4823824 9.647648e-01 5.176176e-01
[7,] 0.5204577 9.590847e-01 4.795423e-01
[8,] 0.9975831 4.833715e-03 2.416857e-03
[9,] 0.9983649 3.270198e-03 1.635099e-03
[10,] 0.9990661 1.867881e-03 9.339403e-04
[11,] 0.9999417 1.166014e-04 5.830071e-05
[12,] 0.9999284 1.432670e-04 7.163351e-05
[13,] 0.9999484 1.032819e-04 5.164097e-05
[14,] 0.9999998 4.352993e-07 2.176497e-07
[15,] 0.9999999 1.593661e-07 7.968307e-08
[16,] 1.0000000 4.254861e-08 2.127430e-08
[17,] 1.0000000 1.491929e-08 7.459646e-09
[18,] 1.0000000 1.490551e-08 7.452756e-09
[19,] 1.0000000 5.641706e-09 2.820853e-09
[20,] 1.0000000 5.107233e-09 2.553617e-09
[21,] 1.0000000 2.421489e-09 1.210744e-09
[22,] 1.0000000 5.735756e-09 2.867878e-09
[23,] 1.0000000 1.050295e-08 5.251475e-09
[24,] 1.0000000 6.418254e-09 3.209127e-09
[25,] 1.0000000 2.361251e-08 1.180626e-08
[26,] 1.0000000 3.865366e-08 1.932683e-08
[27,] 0.9999999 1.116002e-07 5.580012e-08
[28,] 0.9999999 2.384446e-07 1.192223e-07
[29,] 0.9999999 2.767761e-07 1.383881e-07
[30,] 0.9999997 5.959096e-07 2.979548e-07
[31,] 0.9999997 6.243545e-07 3.121772e-07
[32,] 0.9999993 1.423666e-06 7.118331e-07
[33,] 0.9999990 1.925934e-06 9.629669e-07
[34,] 0.9999968 6.365335e-06 3.182668e-06
[35,] 0.9999898 2.044330e-05 1.022165e-05
[36,] 0.9999920 1.600661e-05 8.003307e-06
[37,] 0.9999865 2.706983e-05 1.353492e-05
[38,] 0.9999687 6.269130e-05 3.134565e-05
[39,] 0.9999078 1.843822e-04 9.219108e-05
[40,] 0.9999043 1.913735e-04 9.568675e-05
[41,] 0.9997709 4.581708e-04 2.290854e-04
[42,] 0.9997392 5.216779e-04 2.608390e-04
[43,] 0.9990785 1.843092e-03 9.215462e-04
[44,] 0.9954154 9.169179e-03 4.584589e-03
[45,] 0.9808190 3.836208e-02 1.918104e-02
> postscript(file="/var/wessaorg/rcomp/tmp/1poor1322159742.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/wessaorg/rcomp/tmp/2mwk11322159742.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/wessaorg/rcomp/tmp/3xw7y1322159742.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/wessaorg/rcomp/tmp/48mvj1322159742.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/wessaorg/rcomp/tmp/5pim71322159742.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 = 60
Frequency = 1
1 2 3 4 5 6
-2.174971587 -0.837559832 -0.655700575 -0.119347107 0.106954705 0.822590706
7 8 9 10 11 12
0.915300760 0.893344156 0.561552594 0.569047727 0.605964388 0.709917186
13 14 15 16 17 18
0.770717523 -0.161993275 0.504454215 0.059999562 -0.338378098 0.307402645
19 20 21 22 23 24
-0.091492185 -0.291472911 -0.214828256 -0.711721249 -0.610412968 -0.336769879
25 26 27 28 29 30
-0.382744457 -0.072442304 0.213501280 0.485935945 -0.176896155 0.048692951
31 32 33 34 35 36
0.258914199 0.013670499 0.260330406 0.317713707 0.070905205 -0.037425822
37 38 39 40 41 42
-0.542004294 -0.247572832 -0.519010623 -0.639162405 -0.240493964 0.164467379
43 44 45 46 47 48
0.170104653 0.469321694 0.354753267 0.194029974 -0.474611335 -0.200725287
49 50 51 52 53 54
-0.234881236 -0.247989964 -0.048574260 -0.188574893 0.008380593 0.193001110
55 56 57 58 59 60
0.097105488 -0.076229152 -0.280994444 0.357641256 0.218118772 0.431146805
> postscript(file="/var/wessaorg/rcomp/tmp/66cln1322159742.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -2.174971587 NA
1 -0.837559832 -2.174971587
2 -0.655700575 -0.837559832
3 -0.119347107 -0.655700575
4 0.106954705 -0.119347107
5 0.822590706 0.106954705
6 0.915300760 0.822590706
7 0.893344156 0.915300760
8 0.561552594 0.893344156
9 0.569047727 0.561552594
10 0.605964388 0.569047727
11 0.709917186 0.605964388
12 0.770717523 0.709917186
13 -0.161993275 0.770717523
14 0.504454215 -0.161993275
15 0.059999562 0.504454215
16 -0.338378098 0.059999562
17 0.307402645 -0.338378098
18 -0.091492185 0.307402645
19 -0.291472911 -0.091492185
20 -0.214828256 -0.291472911
21 -0.711721249 -0.214828256
22 -0.610412968 -0.711721249
23 -0.336769879 -0.610412968
24 -0.382744457 -0.336769879
25 -0.072442304 -0.382744457
26 0.213501280 -0.072442304
27 0.485935945 0.213501280
28 -0.176896155 0.485935945
29 0.048692951 -0.176896155
30 0.258914199 0.048692951
31 0.013670499 0.258914199
32 0.260330406 0.013670499
33 0.317713707 0.260330406
34 0.070905205 0.317713707
35 -0.037425822 0.070905205
36 -0.542004294 -0.037425822
37 -0.247572832 -0.542004294
38 -0.519010623 -0.247572832
39 -0.639162405 -0.519010623
40 -0.240493964 -0.639162405
41 0.164467379 -0.240493964
42 0.170104653 0.164467379
43 0.469321694 0.170104653
44 0.354753267 0.469321694
45 0.194029974 0.354753267
46 -0.474611335 0.194029974
47 -0.200725287 -0.474611335
48 -0.234881236 -0.200725287
49 -0.247989964 -0.234881236
50 -0.048574260 -0.247989964
51 -0.188574893 -0.048574260
52 0.008380593 -0.188574893
53 0.193001110 0.008380593
54 0.097105488 0.193001110
55 -0.076229152 0.097105488
56 -0.280994444 -0.076229152
57 0.357641256 -0.280994444
58 0.218118772 0.357641256
59 0.431146805 0.218118772
60 NA 0.431146805
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.837559832 -2.174971587
[2,] -0.655700575 -0.837559832
[3,] -0.119347107 -0.655700575
[4,] 0.106954705 -0.119347107
[5,] 0.822590706 0.106954705
[6,] 0.915300760 0.822590706
[7,] 0.893344156 0.915300760
[8,] 0.561552594 0.893344156
[9,] 0.569047727 0.561552594
[10,] 0.605964388 0.569047727
[11,] 0.709917186 0.605964388
[12,] 0.770717523 0.709917186
[13,] -0.161993275 0.770717523
[14,] 0.504454215 -0.161993275
[15,] 0.059999562 0.504454215
[16,] -0.338378098 0.059999562
[17,] 0.307402645 -0.338378098
[18,] -0.091492185 0.307402645
[19,] -0.291472911 -0.091492185
[20,] -0.214828256 -0.291472911
[21,] -0.711721249 -0.214828256
[22,] -0.610412968 -0.711721249
[23,] -0.336769879 -0.610412968
[24,] -0.382744457 -0.336769879
[25,] -0.072442304 -0.382744457
[26,] 0.213501280 -0.072442304
[27,] 0.485935945 0.213501280
[28,] -0.176896155 0.485935945
[29,] 0.048692951 -0.176896155
[30,] 0.258914199 0.048692951
[31,] 0.013670499 0.258914199
[32,] 0.260330406 0.013670499
[33,] 0.317713707 0.260330406
[34,] 0.070905205 0.317713707
[35,] -0.037425822 0.070905205
[36,] -0.542004294 -0.037425822
[37,] -0.247572832 -0.542004294
[38,] -0.519010623 -0.247572832
[39,] -0.639162405 -0.519010623
[40,] -0.240493964 -0.639162405
[41,] 0.164467379 -0.240493964
[42,] 0.170104653 0.164467379
[43,] 0.469321694 0.170104653
[44,] 0.354753267 0.469321694
[45,] 0.194029974 0.354753267
[46,] -0.474611335 0.194029974
[47,] -0.200725287 -0.474611335
[48,] -0.234881236 -0.200725287
[49,] -0.247989964 -0.234881236
[50,] -0.048574260 -0.247989964
[51,] -0.188574893 -0.048574260
[52,] 0.008380593 -0.188574893
[53,] 0.193001110 0.008380593
[54,] 0.097105488 0.193001110
[55,] -0.076229152 0.097105488
[56,] -0.280994444 -0.076229152
[57,] 0.357641256 -0.280994444
[58,] 0.218118772 0.357641256
[59,] 0.431146805 0.218118772
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.837559832 -2.174971587
2 -0.655700575 -0.837559832
3 -0.119347107 -0.655700575
4 0.106954705 -0.119347107
5 0.822590706 0.106954705
6 0.915300760 0.822590706
7 0.893344156 0.915300760
8 0.561552594 0.893344156
9 0.569047727 0.561552594
10 0.605964388 0.569047727
11 0.709917186 0.605964388
12 0.770717523 0.709917186
13 -0.161993275 0.770717523
14 0.504454215 -0.161993275
15 0.059999562 0.504454215
16 -0.338378098 0.059999562
17 0.307402645 -0.338378098
18 -0.091492185 0.307402645
19 -0.291472911 -0.091492185
20 -0.214828256 -0.291472911
21 -0.711721249 -0.214828256
22 -0.610412968 -0.711721249
23 -0.336769879 -0.610412968
24 -0.382744457 -0.336769879
25 -0.072442304 -0.382744457
26 0.213501280 -0.072442304
27 0.485935945 0.213501280
28 -0.176896155 0.485935945
29 0.048692951 -0.176896155
30 0.258914199 0.048692951
31 0.013670499 0.258914199
32 0.260330406 0.013670499
33 0.317713707 0.260330406
34 0.070905205 0.317713707
35 -0.037425822 0.070905205
36 -0.542004294 -0.037425822
37 -0.247572832 -0.542004294
38 -0.519010623 -0.247572832
39 -0.639162405 -0.519010623
40 -0.240493964 -0.639162405
41 0.164467379 -0.240493964
42 0.170104653 0.164467379
43 0.469321694 0.170104653
44 0.354753267 0.469321694
45 0.194029974 0.354753267
46 -0.474611335 0.194029974
47 -0.200725287 -0.474611335
48 -0.234881236 -0.200725287
49 -0.247989964 -0.234881236
50 -0.048574260 -0.247989964
51 -0.188574893 -0.048574260
52 0.008380593 -0.188574893
53 0.193001110 0.008380593
54 0.097105488 0.193001110
55 -0.076229152 0.097105488
56 -0.280994444 -0.076229152
57 0.357641256 -0.280994444
58 0.218118772 0.357641256
59 0.431146805 0.218118772
> 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/wessaorg/rcomp/tmp/7cbbw1322159742.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/wessaorg/rcomp/tmp/86dnv1322159742.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/wessaorg/rcomp/tmp/92bno1322159742.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/wessaorg/rcomp/tmp/1087p51322159742.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/114x3c1322159742.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/wessaorg/rcomp/tmp/12kzw91322159742.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/wessaorg/rcomp/tmp/13ijhp1322159742.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/wessaorg/rcomp/tmp/14v3c71322159742.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/wessaorg/rcomp/tmp/15en3k1322159742.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/wessaorg/rcomp/tmp/161qg31322159742.tab")
+ }
>
> try(system("convert tmp/1poor1322159742.ps tmp/1poor1322159742.png",intern=TRUE))
character(0)
> try(system("convert tmp/2mwk11322159742.ps tmp/2mwk11322159742.png",intern=TRUE))
character(0)
> try(system("convert tmp/3xw7y1322159742.ps tmp/3xw7y1322159742.png",intern=TRUE))
character(0)
> try(system("convert tmp/48mvj1322159742.ps tmp/48mvj1322159742.png",intern=TRUE))
character(0)
> try(system("convert tmp/5pim71322159742.ps tmp/5pim71322159742.png",intern=TRUE))
character(0)
> try(system("convert tmp/66cln1322159742.ps tmp/66cln1322159742.png",intern=TRUE))
character(0)
> try(system("convert tmp/7cbbw1322159742.ps tmp/7cbbw1322159742.png",intern=TRUE))
character(0)
> try(system("convert tmp/86dnv1322159742.ps tmp/86dnv1322159742.png",intern=TRUE))
character(0)
> try(system("convert tmp/92bno1322159742.ps tmp/92bno1322159742.png",intern=TRUE))
character(0)
> try(system("convert tmp/1087p51322159742.ps tmp/1087p51322159742.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.308 0.491 3.811