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 = 'No 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
1 63.46 60635408600 7175.79 8450000
2 64.98 64463389370 7533.41 8557000
3 65.48 66750859382 7749.11 8614000
4 66.35 69032584256 7990.47 8639369
5 66.80 72841303907 8393.42 8678386
6 68.00 72838685608 8343.11 8730405
7 68.37 75887977454 8645.37 8777873
8 68.63 78936158550 8950.31 8819380
9 68.58 81986612544 9244.73 8868475
10 68.87 85038924238 9529.40 8923845
11 69.24 87328862790 9714.96 8989111
12 69.93 85794984681 9477.27 9052707
13 70.33 88082838620 9675.47 9103729
14 69.65 91885510657 10076.60 9118700
15 70.52 96355573324 10512.51 9165800
16 70.25 101321342609 10991.21 9218400
17 70.06 105791395837 11396.13 9283100
18 70.73 113241494103 12089.41 9367000
19 70.58 117215849653 12406.29 9448100
20 70.65 120940898881 12720.18 9507800
21 70.94 125658810316 13149.04 9556500
22 70.63 130873880201 13647.20 9589800
23 70.71 139583536623 14520.74 9612700
24 70.97 148226520849 15379.71 9637800
25 71.10 153789461350 15899.66 9672500
26 71.44 161871513310 16672.14 9709100
27 71.65 171781295610 17639.58 9738400
28 72.01 178996585758 18325.17 9767800
29 72.00 176620960002 18032.12 9794800
30 72.15 186604680395 19019.94 9811000
31 72.80 187772917033 19117.97 9821800
32 72.75 193109744879 19645.54 9829700
33 73.24 197630528464 20090.12 9837200
34 73.29 206483683756 20969.62 9846800
35 73.70 203906590622 20696.13 9852400
36 73.93 206783719069 20979.85 9856303
37 73.93 206759082202 20979.01 9855520
38 74.43 211878483671 21498.94 9855300
39 74.54 214009149959 21708.75 9858200
40 74.74 217203709135 22024.75 9861800
41 75.35 222331811923 22525.56 9870200
42 75.66 232825724880 23555.82 9884000
43 75.73 241180274039 24269.23 9937697
44 76.14 248494123822 24925.91 9969310
45 76.30 253049202253 25293.57 10004487
46 76.46 256922518700 25575.57 10045622
47 76.47 254451243639 25229.60 10085426
48 76.82 262662316801 25947.30 10122914
49 76.97 268926171540 26480.95 10155459
50 77.31 272037080259 26725.50 10178934
51 77.53 281118338865 27561.20 10199787
52 77.62 286389613939 28030.61 10217030
53 77.80 296190694319 28937.15 10235655
54 77.91 307294826962 29940.20 10263618
55 78.22 309697986636 30092.08 10291679
56 78.32 314369521231 30485.88 10311970
57 78.47 317533640477 30736.53 10330824
58 79.12 327005697521 31600.02 10348276
59 79.21 332458473876 32077.00 10364388
60 79.55 339794119726 32738.41 10379067
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) GDP Inkomen Populatie
4.326e+01 3.655e-11 -6.945e-05 2.509e-06
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.7183 -0.5358 0.1204 0.4635 1.6817
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.326e+01 5.540e+00 7.807 1.6e-10 ***
GDP 3.655e-11 6.163e-11 0.593 0.555479
Inkomen -6.945e-05 6.769e-04 -0.103 0.918642
Populatie 2.509e-06 6.710e-07 3.740 0.000435 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8335 on 56 degrees of freedom
Multiple R-squared: 0.9578, Adjusted R-squared: 0.9556
F-statistic: 423.9 on 3 and 56 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.2106343 4.212686e-01 7.893657e-01
[2,] 0.2113188 4.226376e-01 7.886812e-01
[3,] 0.3358066 6.716132e-01 6.641934e-01
[4,] 0.3334964 6.669928e-01 6.665036e-01
[5,] 0.3046183 6.092367e-01 6.953817e-01
[6,] 0.2097174 4.194349e-01 7.902826e-01
[7,] 0.1595751 3.191503e-01 8.404249e-01
[8,] 0.1163143 2.326287e-01 8.836857e-01
[9,] 0.3759185 7.518370e-01 6.240815e-01
[10,] 0.3239842 6.479684e-01 6.760158e-01
[11,] 0.2608690 5.217380e-01 7.391310e-01
[12,] 0.4445311 8.890622e-01 5.554689e-01
[13,] 0.4361025 8.722050e-01 5.638975e-01
[14,] 0.4912038 9.824076e-01 5.087962e-01
[15,] 0.8079611 3.840778e-01 1.920389e-01
[16,] 0.9197322 1.605356e-01 8.026781e-02
[17,] 0.9782065 4.358697e-02 2.179349e-02
[18,] 0.9954648 9.070349e-03 4.535174e-03
[19,] 0.9989968 2.006400e-03 1.003200e-03
[20,] 0.9998660 2.679904e-04 1.339952e-04
[21,] 0.9999081 1.837550e-04 9.187748e-05
[22,] 0.9999451 1.097649e-04 5.488245e-05
[23,] 0.9999747 5.063039e-05 2.531520e-05
[24,] 0.9999894 2.127852e-05 1.063926e-05
[25,] 0.9999956 8.882773e-06 4.441387e-06
[26,] 0.9999965 7.043213e-06 3.521606e-06
[27,] 0.9999965 6.936464e-06 3.468232e-06
[28,] 0.9999999 2.537171e-07 1.268586e-07
[29,] 0.9999999 1.840189e-07 9.200943e-08
[30,] 0.9999999 1.601336e-07 8.006679e-08
[31,] 1.0000000 6.961932e-08 3.480966e-08
[32,] 1.0000000 9.942761e-08 4.971380e-08
[33,] 1.0000000 7.586013e-08 3.793007e-08
[34,] 1.0000000 6.504329e-09 3.252165e-09
[35,] 1.0000000 1.141373e-08 5.706864e-09
[36,] 1.0000000 5.343279e-08 2.671639e-08
[37,] 1.0000000 3.560088e-08 1.780044e-08
[38,] 1.0000000 8.210121e-08 4.105061e-08
[39,] 0.9999999 2.174508e-07 1.087254e-07
[40,] 0.9999996 7.507763e-07 3.753882e-07
[41,] 0.9999991 1.833170e-06 9.165851e-07
[42,] 0.9999967 6.672052e-06 3.336026e-06
[43,] 0.9999981 3.841014e-06 1.920507e-06
[44,] 0.9999914 1.713746e-05 8.568732e-06
[45,] 0.9999227 1.546155e-04 7.730775e-05
[46,] 0.9993383 1.323488e-03 6.617440e-04
[47,] 0.9976031 4.793775e-03 2.396887e-03
> postscript(file="/var/wessaorg/rcomp/tmp/1cy2x1322154598.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/2vysr1322154598.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/3260y1322154598.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/4ir2a1322154598.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/5vtwq1322154598.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.718255757 -1.581836664 -1.293499229 -0.553797674 -0.312935878 0.753133894
7 8 9 10 11 12
0.913555348 0.979161526 0.714913339 0.774174165 0.909586094 1.479563007
13 14 15 16 17 18
1.681671619 0.852966901 1.471661301 0.921407027 0.433785575 0.669084462
19 20 21 22 23 24
0.192314725 -0.001850666 0.023279662 -0.526305988 -0.761459773 -0.820708929
25 26 27 28 29 30
-0.945009813 -0.938619492 -1.097178951 -1.027073981 -1.038343246 -1.225317933
31 32 33 34 35 36
-0.638312086 -0.866569244 -0.529758045 -0.766369347 -0.295216729 -0.160471884
37 38 39 40 41 42
-0.157664894 0.191870420 0.231283963 0.327428153 0.763687218 0.727033587
43 44 45 46 47 48
0.406459142 0.515400495 0.446165625 0.380951115 0.357373417 0.363014739
49 50 51 52 53 54
0.239452653 0.423819524 0.317591126 0.204246207 0.042217240 -0.254170709
55 56 57 58 59 60
-0.091877680 -0.186200248 -0.181759181 0.138191973 0.021576491 0.102542287
> postscript(file="/var/wessaorg/rcomp/tmp/6v70h1322154598.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.718255757 NA
1 -1.581836664 -2.718255757
2 -1.293499229 -1.581836664
3 -0.553797674 -1.293499229
4 -0.312935878 -0.553797674
5 0.753133894 -0.312935878
6 0.913555348 0.753133894
7 0.979161526 0.913555348
8 0.714913339 0.979161526
9 0.774174165 0.714913339
10 0.909586094 0.774174165
11 1.479563007 0.909586094
12 1.681671619 1.479563007
13 0.852966901 1.681671619
14 1.471661301 0.852966901
15 0.921407027 1.471661301
16 0.433785575 0.921407027
17 0.669084462 0.433785575
18 0.192314725 0.669084462
19 -0.001850666 0.192314725
20 0.023279662 -0.001850666
21 -0.526305988 0.023279662
22 -0.761459773 -0.526305988
23 -0.820708929 -0.761459773
24 -0.945009813 -0.820708929
25 -0.938619492 -0.945009813
26 -1.097178951 -0.938619492
27 -1.027073981 -1.097178951
28 -1.038343246 -1.027073981
29 -1.225317933 -1.038343246
30 -0.638312086 -1.225317933
31 -0.866569244 -0.638312086
32 -0.529758045 -0.866569244
33 -0.766369347 -0.529758045
34 -0.295216729 -0.766369347
35 -0.160471884 -0.295216729
36 -0.157664894 -0.160471884
37 0.191870420 -0.157664894
38 0.231283963 0.191870420
39 0.327428153 0.231283963
40 0.763687218 0.327428153
41 0.727033587 0.763687218
42 0.406459142 0.727033587
43 0.515400495 0.406459142
44 0.446165625 0.515400495
45 0.380951115 0.446165625
46 0.357373417 0.380951115
47 0.363014739 0.357373417
48 0.239452653 0.363014739
49 0.423819524 0.239452653
50 0.317591126 0.423819524
51 0.204246207 0.317591126
52 0.042217240 0.204246207
53 -0.254170709 0.042217240
54 -0.091877680 -0.254170709
55 -0.186200248 -0.091877680
56 -0.181759181 -0.186200248
57 0.138191973 -0.181759181
58 0.021576491 0.138191973
59 0.102542287 0.021576491
60 NA 0.102542287
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.581836664 -2.718255757
[2,] -1.293499229 -1.581836664
[3,] -0.553797674 -1.293499229
[4,] -0.312935878 -0.553797674
[5,] 0.753133894 -0.312935878
[6,] 0.913555348 0.753133894
[7,] 0.979161526 0.913555348
[8,] 0.714913339 0.979161526
[9,] 0.774174165 0.714913339
[10,] 0.909586094 0.774174165
[11,] 1.479563007 0.909586094
[12,] 1.681671619 1.479563007
[13,] 0.852966901 1.681671619
[14,] 1.471661301 0.852966901
[15,] 0.921407027 1.471661301
[16,] 0.433785575 0.921407027
[17,] 0.669084462 0.433785575
[18,] 0.192314725 0.669084462
[19,] -0.001850666 0.192314725
[20,] 0.023279662 -0.001850666
[21,] -0.526305988 0.023279662
[22,] -0.761459773 -0.526305988
[23,] -0.820708929 -0.761459773
[24,] -0.945009813 -0.820708929
[25,] -0.938619492 -0.945009813
[26,] -1.097178951 -0.938619492
[27,] -1.027073981 -1.097178951
[28,] -1.038343246 -1.027073981
[29,] -1.225317933 -1.038343246
[30,] -0.638312086 -1.225317933
[31,] -0.866569244 -0.638312086
[32,] -0.529758045 -0.866569244
[33,] -0.766369347 -0.529758045
[34,] -0.295216729 -0.766369347
[35,] -0.160471884 -0.295216729
[36,] -0.157664894 -0.160471884
[37,] 0.191870420 -0.157664894
[38,] 0.231283963 0.191870420
[39,] 0.327428153 0.231283963
[40,] 0.763687218 0.327428153
[41,] 0.727033587 0.763687218
[42,] 0.406459142 0.727033587
[43,] 0.515400495 0.406459142
[44,] 0.446165625 0.515400495
[45,] 0.380951115 0.446165625
[46,] 0.357373417 0.380951115
[47,] 0.363014739 0.357373417
[48,] 0.239452653 0.363014739
[49,] 0.423819524 0.239452653
[50,] 0.317591126 0.423819524
[51,] 0.204246207 0.317591126
[52,] 0.042217240 0.204246207
[53,] -0.254170709 0.042217240
[54,] -0.091877680 -0.254170709
[55,] -0.186200248 -0.091877680
[56,] -0.181759181 -0.186200248
[57,] 0.138191973 -0.181759181
[58,] 0.021576491 0.138191973
[59,] 0.102542287 0.021576491
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.581836664 -2.718255757
2 -1.293499229 -1.581836664
3 -0.553797674 -1.293499229
4 -0.312935878 -0.553797674
5 0.753133894 -0.312935878
6 0.913555348 0.753133894
7 0.979161526 0.913555348
8 0.714913339 0.979161526
9 0.774174165 0.714913339
10 0.909586094 0.774174165
11 1.479563007 0.909586094
12 1.681671619 1.479563007
13 0.852966901 1.681671619
14 1.471661301 0.852966901
15 0.921407027 1.471661301
16 0.433785575 0.921407027
17 0.669084462 0.433785575
18 0.192314725 0.669084462
19 -0.001850666 0.192314725
20 0.023279662 -0.001850666
21 -0.526305988 0.023279662
22 -0.761459773 -0.526305988
23 -0.820708929 -0.761459773
24 -0.945009813 -0.820708929
25 -0.938619492 -0.945009813
26 -1.097178951 -0.938619492
27 -1.027073981 -1.097178951
28 -1.038343246 -1.027073981
29 -1.225317933 -1.038343246
30 -0.638312086 -1.225317933
31 -0.866569244 -0.638312086
32 -0.529758045 -0.866569244
33 -0.766369347 -0.529758045
34 -0.295216729 -0.766369347
35 -0.160471884 -0.295216729
36 -0.157664894 -0.160471884
37 0.191870420 -0.157664894
38 0.231283963 0.191870420
39 0.327428153 0.231283963
40 0.763687218 0.327428153
41 0.727033587 0.763687218
42 0.406459142 0.727033587
43 0.515400495 0.406459142
44 0.446165625 0.515400495
45 0.380951115 0.446165625
46 0.357373417 0.380951115
47 0.363014739 0.357373417
48 0.239452653 0.363014739
49 0.423819524 0.239452653
50 0.317591126 0.423819524
51 0.204246207 0.317591126
52 0.042217240 0.204246207
53 -0.254170709 0.042217240
54 -0.091877680 -0.254170709
55 -0.186200248 -0.091877680
56 -0.181759181 -0.186200248
57 0.138191973 -0.181759181
58 0.021576491 0.138191973
59 0.102542287 0.021576491
> 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/7wu3n1322154598.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/8m7cz1322154598.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/95fw01322154598.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/10dprf1322154598.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/11vsed1322154598.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/12kgog1322154598.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/139k651322154598.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/14s6no1322154598.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/151rcn1322154598.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/16ypn41322154598.tab")
+ }
>
> try(system("convert tmp/1cy2x1322154598.ps tmp/1cy2x1322154598.png",intern=TRUE))
character(0)
> try(system("convert tmp/2vysr1322154598.ps tmp/2vysr1322154598.png",intern=TRUE))
character(0)
> try(system("convert tmp/3260y1322154598.ps tmp/3260y1322154598.png",intern=TRUE))
character(0)
> try(system("convert tmp/4ir2a1322154598.ps tmp/4ir2a1322154598.png",intern=TRUE))
character(0)
> try(system("convert tmp/5vtwq1322154598.ps tmp/5vtwq1322154598.png",intern=TRUE))
character(0)
> try(system("convert tmp/6v70h1322154598.ps tmp/6v70h1322154598.png",intern=TRUE))
character(0)
> try(system("convert tmp/7wu3n1322154598.ps tmp/7wu3n1322154598.png",intern=TRUE))
character(0)
> try(system("convert tmp/8m7cz1322154598.ps tmp/8m7cz1322154598.png",intern=TRUE))
character(0)
> try(system("convert tmp/95fw01322154598.ps tmp/95fw01322154598.png",intern=TRUE))
character(0)
> try(system("convert tmp/10dprf1322154598.ps tmp/10dprf1322154598.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.330 0.540 3.984