R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-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(132838
+ ,312991
+ ,5599
+ ,47645
+ ,15545
+ ,575093
+ ,129842
+ ,301647
+ ,5234
+ ,45970
+ ,15001
+ ,557560
+ ,129694
+ ,305353
+ ,5279
+ ,48069
+ ,14961
+ ,564478
+ ,130080
+ ,313665
+ ,5391
+ ,53080
+ ,15245
+ ,580523
+ ,131496
+ ,322402
+ ,5280
+ ,57896
+ ,15656
+ ,596594
+ ,131556
+ ,318280
+ ,5173
+ ,54344
+ ,15577
+ ,586570
+ ,128925
+ ,292852
+ ,4724
+ ,40482
+ ,14630
+ ,536214
+ ,127836
+ ,287481
+ ,4554
+ ,37110
+ ,14336
+ ,523597
+ ,129164
+ ,295210
+ ,4713
+ ,39263
+ ,14834
+ ,536535
+ ,129531
+ ,295650
+ ,4811
+ ,38889
+ ,14921
+ ,536322
+ ,128548
+ ,292919
+ ,4668
+ ,39593
+ ,14707
+ ,532638
+ ,127330
+ ,290649
+ ,4516
+ ,39305
+ ,14516
+ ,528222
+ ,123815
+ ,281687
+ ,4203
+ ,40560
+ ,14055
+ ,516141
+ ,124393
+ ,270336
+ ,4016
+ ,38306
+ ,13493
+ ,501866
+ ,123707
+ ,271420
+ ,3993
+ ,40911
+ ,13528
+ ,506174
+ ,123736
+ ,278183
+ ,3971
+ ,44700
+ ,13719
+ ,517945
+ ,124507
+ ,284913
+ ,3838
+ ,50328
+ ,14170
+ ,533590
+ ,125005
+ ,283487
+ ,3891
+ ,47499
+ ,14009
+ ,528379
+ ,121383
+ ,256677
+ ,3306
+ ,34446
+ ,13159
+ ,477580
+ ,121200
+ ,252945
+ ,3235
+ ,31434
+ ,12927
+ ,469357
+ ,125249
+ ,264963
+ ,3404
+ ,34066
+ ,13510
+ ,490243
+ ,125253
+ ,265988
+ ,3400
+ ,35044
+ ,13520
+ ,492622
+ ,127977
+ ,274857
+ ,3447
+ ,37040
+ ,14089
+ ,507561
+ ,128984
+ ,279650
+ ,3431
+ ,38706
+ ,14251
+ ,516922
+ ,126770
+ ,276715
+ ,3321
+ ,40430
+ ,13980
+ ,514258
+ ,126448
+ ,273887
+ ,3189
+ ,39613
+ ,13715
+ ,509846
+ ,127845
+ ,282308
+ ,3256
+ ,44236
+ ,14112
+ ,527070
+ ,128818
+ ,289847
+ ,3290
+ ,47859
+ ,14289
+ ,541657
+ ,132127
+ ,301101
+ ,3475
+ ,53711
+ ,15020
+ ,564591
+ ,132338
+ ,297008
+ ,3454
+ ,50352
+ ,14860
+ ,555362
+ ,126645
+ ,268909
+ ,2806
+ ,36142
+ ,13800
+ ,498662
+ ,130625
+ ,278383
+ ,2777
+ ,34819
+ ,14431
+ ,511038
+ ,133506
+ ,286226
+ ,2865
+ ,37353
+ ,14944
+ ,525919
+ ,135277
+ ,288936
+ ,2924
+ ,37550
+ ,15083
+ ,531673
+ ,137664
+ ,298953
+ ,3011
+ ,40462
+ ,15707
+ ,548854
+ ,139821
+ ,305837
+ ,3099
+ ,41753
+ ,15954
+ ,560576
+ ,138440
+ ,301979
+ ,2988
+ ,43437
+ ,15631
+ ,557274
+ ,139879
+ ,306281
+ ,3032
+ ,44784
+ ,15813
+ ,565742
+ ,142256
+ ,317057
+ ,3131
+ ,49537
+ ,16356
+ ,587625
+ ,146322
+ ,334780
+ ,3343
+ ,54974
+ ,17086
+ ,619916
+ ,146389
+ ,335895
+ ,3275
+ ,58535
+ ,17302
+ ,625809
+ ,147841
+ ,333874
+ ,3243
+ ,54762
+ ,17247
+ ,619567
+ ,146449
+ ,311028
+ ,2897
+ ,40738
+ ,16398
+ ,572942
+ ,147960
+ ,311767
+ ,2818
+ ,38052
+ ,16590
+ ,572775
+ ,148487
+ ,312575
+ ,2836
+ ,38436
+ ,16673
+ ,574205
+ ,149802
+ ,315040
+ ,2721
+ ,36993
+ ,16962
+ ,579799
+ ,151387
+ ,320325
+ ,2742
+ ,39056
+ ,17278
+ ,590072
+ ,151936
+ ,321178
+ ,2707
+ ,39996
+ ,17224
+ ,593408)
+ ,dim=c(6
+ ,48)
+ ,dimnames=list(c('Basisonderwijs(lager_1ste_graad_secundair)'
+ ,'Secundair_onderwijs(2de
+ ,3de
+ ,4de_graad)'
+ ,'Duaal_onderwijs'
+ ,'Hoger_onderwijs(Bachelor)'
+ ,'Leercontract'
+ ,'Werkloosheid_totaal')
+ ,1:48))
> y <- array(NA,dim=c(6,48),dimnames=list(c('Basisonderwijs(lager_1ste_graad_secundair)','Secundair_onderwijs(2de,3de,4de_graad)','Duaal_onderwijs','Hoger_onderwijs(Bachelor)','Leercontract','Werkloosheid_totaal'),1:48))
> 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 = '6'
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, as.Date.numeric
> 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
Werkloosheid_totaal Basisonderwijs(lager_1ste_graad_secundair)
1 575093 132838
2 557560 129842
3 564478 129694
4 580523 130080
5 596594 131496
6 586570 131556
7 536214 128925
8 523597 127836
9 536535 129164
10 536322 129531
11 532638 128548
12 528222 127330
13 516141 123815
14 501866 124393
15 506174 123707
16 517945 123736
17 533590 124507
18 528379 125005
19 477580 121383
20 469357 121200
21 490243 125249
22 492622 125253
23 507561 127977
24 516922 128984
25 514258 126770
26 509846 126448
27 527070 127845
28 541657 128818
29 564591 132127
30 555362 132338
31 498662 126645
32 511038 130625
33 525919 133506
34 531673 135277
35 548854 137664
36 560576 139821
37 557274 138440
38 565742 139879
39 587625 142256
40 619916 146322
41 625809 146389
42 619567 147841
43 572942 146449
44 572775 147960
45 574205 148487
46 579799 149802
47 590072 151387
48 593408 151936
Secundair_onderwijs(2de,3de,4de_graad) Duaal_onderwijs
1 312991 5599
2 301647 5234
3 305353 5279
4 313665 5391
5 322402 5280
6 318280 5173
7 292852 4724
8 287481 4554
9 295210 4713
10 295650 4811
11 292919 4668
12 290649 4516
13 281687 4203
14 270336 4016
15 271420 3993
16 278183 3971
17 284913 3838
18 283487 3891
19 256677 3306
20 252945 3235
21 264963 3404
22 265988 3400
23 274857 3447
24 279650 3431
25 276715 3321
26 273887 3189
27 282308 3256
28 289847 3290
29 301101 3475
30 297008 3454
31 268909 2806
32 278383 2777
33 286226 2865
34 288936 2924
35 298953 3011
36 305837 3099
37 301979 2988
38 306281 3032
39 317057 3131
40 334780 3343
41 335895 3275
42 333874 3243
43 311028 2897
44 311767 2818
45 312575 2836
46 315040 2721
47 320325 2742
48 321178 2707
Hoger_onderwijs(Bachelor) Leercontract
1 47645 15545
2 45970 15001
3 48069 14961
4 53080 15245
5 57896 15656
6 54344 15577
7 40482 14630
8 37110 14336
9 39263 14834
10 38889 14921
11 39593 14707
12 39305 14516
13 40560 14055
14 38306 13493
15 40911 13528
16 44700 13719
17 50328 14170
18 47499 14009
19 34446 13159
20 31434 12927
21 34066 13510
22 35044 13520
23 37040 14089
24 38706 14251
25 40430 13980
26 39613 13715
27 44236 14112
28 47859 14289
29 53711 15020
30 50352 14860
31 36142 13800
32 34819 14431
33 37353 14944
34 37550 15083
35 40462 15707
36 41753 15954
37 43437 15631
38 44784 15813
39 49537 16356
40 54974 17086
41 58535 17302
42 54762 17247
43 40738 16398
44 38052 16590
45 38436 16673
46 36993 16962
47 39056 17278
48 39996 17224
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept)
-22725.320
`Basisonderwijs(lager_1ste_graad_secundair)`
1.722
`Secundair_onderwijs(2de,3de,4de_graad)`
0.997
Duaal_onderwijs
3.394
`Hoger_onderwijs(Bachelor)`
1.500
Leercontract
-2.159
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2127.2 -956.0 -239.6 932.8 2391.2
Coefficients:
Estimate Std. Error t value
(Intercept) -2.273e+04 8.210e+03 -2.768
`Basisonderwijs(lager_1ste_graad_secundair)` 1.722e+00 1.797e-01 9.581
`Secundair_onderwijs(2de,3de,4de_graad)` 9.970e-01 1.226e-01 8.131
Duaal_onderwijs 3.394e+00 7.340e-01 4.624
`Hoger_onderwijs(Bachelor)` 1.500e+00 9.848e-02 15.237
Leercontract -2.159e+00 2.065e+00 -1.046
Pr(>|t|)
(Intercept) 0.00836 **
`Basisonderwijs(lager_1ste_graad_secundair)` 3.96e-12 ***
`Secundair_onderwijs(2de,3de,4de_graad)` 3.67e-10 ***
Duaal_onderwijs 3.57e-05 ***
`Hoger_onderwijs(Bachelor)` < 2e-16 ***
Leercontract 0.30166
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1175 on 42 degrees of freedom
Multiple R-squared: 0.9991, Adjusted R-squared: 0.999
F-statistic: 9694 on 5 and 42 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.016891112 0.0337822249 0.9831088876
[2,] 0.031840351 0.0636807028 0.9681596486
[3,] 0.075158415 0.1503168305 0.9248415848
[4,] 0.038722857 0.0774457142 0.9612771429
[5,] 0.019042316 0.0380846316 0.9809576842
[6,] 0.018882139 0.0377642785 0.9811178608
[7,] 0.008356227 0.0167124539 0.9916437731
[8,] 0.003449856 0.0068997129 0.9965501436
[9,] 0.006033269 0.0120665376 0.9939667312
[10,] 0.002532931 0.0050658625 0.9974670688
[11,] 0.028407620 0.0568152397 0.9715923802
[12,] 0.088472663 0.1769453255 0.9115273372
[13,] 0.091966001 0.1839320029 0.9080339986
[14,] 0.065017892 0.1300357842 0.9349821079
[15,] 0.039716330 0.0794326597 0.9602836702
[16,] 0.033335344 0.0666706870 0.9666646565
[17,] 0.069165065 0.1383301310 0.9308349345
[18,] 0.088236384 0.1764727687 0.9117636156
[19,] 0.107971055 0.2159421098 0.8920289451
[20,] 0.163958667 0.3279173332 0.8360413334
[21,] 0.192171818 0.3843436359 0.8078281821
[22,] 0.824195788 0.3516084233 0.1758042117
[23,] 0.987985238 0.0240295243 0.0120147622
[24,] 0.986902839 0.0261943212 0.0130971606
[25,] 0.987662293 0.0246754142 0.0123377071
[26,] 0.999597196 0.0008056087 0.0004028043
[27,] 0.998750678 0.0024986437 0.0012493218
[28,] 0.995724859 0.0085502827 0.0042751414
[29,] 0.995286949 0.0094261013 0.0047130507
[30,] 0.993614548 0.0127709037 0.0063854519
[31,] 0.971259741 0.0574805171 0.0287402585
> postscript(file="/var/wessaorg/rcomp/tmp/1cskb1353353270.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/2fghp1353353270.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/3sjql1353353270.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/4wgqe1353353270.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/5tiui1353353270.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 = 48
Frequency = 1
1 2 3 4 5
124.8620235 1637.6604769 1727.1124212 1534.9301716 495.3626423
6 7 8 9 10
-0.1891302 -196.5112662 -581.9793850 -331.0216955 -1198.2178869
11 12 13 14 15
-1499.9416873 -1019.9987222 70.0121822 -1080.2195842 -426.8055901
16 17 18 19 20
-646.4040972 -57.7786005 -987.3374908 914.7077830 986.9104811
21 22 23 24 25
-344.3056473 -426.3632268 -945.5315003 -192.5329957 1083.1047157
26 27 28 29 30
1146.6922777 1263.0323865 1489.1568312 -324.3933767 -1070.2613160
31 32 33 34 35
1278.1301029 802.0417678 -89.8764086 -282.7074167 -515.5876176
36 37 38 39 40
-1073.1536856 -998.4279162 -1074.5627530 -322.8913735 -3.4485726
41 42 43 44 45
16.7143046 -1059.3955538 -2127.1692536 -919.6256041 -1660.6179372
46 47 48
2391.2083114 2181.7629968 2313.8554150
> postscript(file="/var/wessaorg/rcomp/tmp/6b2zo1353353270.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 = 48
Frequency = 1
lag(myerror, k = 1) myerror
0 124.8620235 NA
1 1637.6604769 124.8620235
2 1727.1124212 1637.6604769
3 1534.9301716 1727.1124212
4 495.3626423 1534.9301716
5 -0.1891302 495.3626423
6 -196.5112662 -0.1891302
7 -581.9793850 -196.5112662
8 -331.0216955 -581.9793850
9 -1198.2178869 -331.0216955
10 -1499.9416873 -1198.2178869
11 -1019.9987222 -1499.9416873
12 70.0121822 -1019.9987222
13 -1080.2195842 70.0121822
14 -426.8055901 -1080.2195842
15 -646.4040972 -426.8055901
16 -57.7786005 -646.4040972
17 -987.3374908 -57.7786005
18 914.7077830 -987.3374908
19 986.9104811 914.7077830
20 -344.3056473 986.9104811
21 -426.3632268 -344.3056473
22 -945.5315003 -426.3632268
23 -192.5329957 -945.5315003
24 1083.1047157 -192.5329957
25 1146.6922777 1083.1047157
26 1263.0323865 1146.6922777
27 1489.1568312 1263.0323865
28 -324.3933767 1489.1568312
29 -1070.2613160 -324.3933767
30 1278.1301029 -1070.2613160
31 802.0417678 1278.1301029
32 -89.8764086 802.0417678
33 -282.7074167 -89.8764086
34 -515.5876176 -282.7074167
35 -1073.1536856 -515.5876176
36 -998.4279162 -1073.1536856
37 -1074.5627530 -998.4279162
38 -322.8913735 -1074.5627530
39 -3.4485726 -322.8913735
40 16.7143046 -3.4485726
41 -1059.3955538 16.7143046
42 -2127.1692536 -1059.3955538
43 -919.6256041 -2127.1692536
44 -1660.6179372 -919.6256041
45 2391.2083114 -1660.6179372
46 2181.7629968 2391.2083114
47 2313.8554150 2181.7629968
48 NA 2313.8554150
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1637.6604769 124.8620235
[2,] 1727.1124212 1637.6604769
[3,] 1534.9301716 1727.1124212
[4,] 495.3626423 1534.9301716
[5,] -0.1891302 495.3626423
[6,] -196.5112662 -0.1891302
[7,] -581.9793850 -196.5112662
[8,] -331.0216955 -581.9793850
[9,] -1198.2178869 -331.0216955
[10,] -1499.9416873 -1198.2178869
[11,] -1019.9987222 -1499.9416873
[12,] 70.0121822 -1019.9987222
[13,] -1080.2195842 70.0121822
[14,] -426.8055901 -1080.2195842
[15,] -646.4040972 -426.8055901
[16,] -57.7786005 -646.4040972
[17,] -987.3374908 -57.7786005
[18,] 914.7077830 -987.3374908
[19,] 986.9104811 914.7077830
[20,] -344.3056473 986.9104811
[21,] -426.3632268 -344.3056473
[22,] -945.5315003 -426.3632268
[23,] -192.5329957 -945.5315003
[24,] 1083.1047157 -192.5329957
[25,] 1146.6922777 1083.1047157
[26,] 1263.0323865 1146.6922777
[27,] 1489.1568312 1263.0323865
[28,] -324.3933767 1489.1568312
[29,] -1070.2613160 -324.3933767
[30,] 1278.1301029 -1070.2613160
[31,] 802.0417678 1278.1301029
[32,] -89.8764086 802.0417678
[33,] -282.7074167 -89.8764086
[34,] -515.5876176 -282.7074167
[35,] -1073.1536856 -515.5876176
[36,] -998.4279162 -1073.1536856
[37,] -1074.5627530 -998.4279162
[38,] -322.8913735 -1074.5627530
[39,] -3.4485726 -322.8913735
[40,] 16.7143046 -3.4485726
[41,] -1059.3955538 16.7143046
[42,] -2127.1692536 -1059.3955538
[43,] -919.6256041 -2127.1692536
[44,] -1660.6179372 -919.6256041
[45,] 2391.2083114 -1660.6179372
[46,] 2181.7629968 2391.2083114
[47,] 2313.8554150 2181.7629968
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1637.6604769 124.8620235
2 1727.1124212 1637.6604769
3 1534.9301716 1727.1124212
4 495.3626423 1534.9301716
5 -0.1891302 495.3626423
6 -196.5112662 -0.1891302
7 -581.9793850 -196.5112662
8 -331.0216955 -581.9793850
9 -1198.2178869 -331.0216955
10 -1499.9416873 -1198.2178869
11 -1019.9987222 -1499.9416873
12 70.0121822 -1019.9987222
13 -1080.2195842 70.0121822
14 -426.8055901 -1080.2195842
15 -646.4040972 -426.8055901
16 -57.7786005 -646.4040972
17 -987.3374908 -57.7786005
18 914.7077830 -987.3374908
19 986.9104811 914.7077830
20 -344.3056473 986.9104811
21 -426.3632268 -344.3056473
22 -945.5315003 -426.3632268
23 -192.5329957 -945.5315003
24 1083.1047157 -192.5329957
25 1146.6922777 1083.1047157
26 1263.0323865 1146.6922777
27 1489.1568312 1263.0323865
28 -324.3933767 1489.1568312
29 -1070.2613160 -324.3933767
30 1278.1301029 -1070.2613160
31 802.0417678 1278.1301029
32 -89.8764086 802.0417678
33 -282.7074167 -89.8764086
34 -515.5876176 -282.7074167
35 -1073.1536856 -515.5876176
36 -998.4279162 -1073.1536856
37 -1074.5627530 -998.4279162
38 -322.8913735 -1074.5627530
39 -3.4485726 -322.8913735
40 16.7143046 -3.4485726
41 -1059.3955538 16.7143046
42 -2127.1692536 -1059.3955538
43 -919.6256041 -2127.1692536
44 -1660.6179372 -919.6256041
45 2391.2083114 -1660.6179372
46 2181.7629968 2391.2083114
47 2313.8554150 2181.7629968
> 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/7qeol1353353270.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/8527g1353353270.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/9tme81353353270.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/10bc9e1353353270.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/11sfri1353353270.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/12z7ge1353353270.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/13q07e1353353270.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/14q2hu1353353270.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/151s8u1353353270.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/166elx1353353270.tab")
+ }
>
> try(system("convert tmp/1cskb1353353270.ps tmp/1cskb1353353270.png",intern=TRUE))
character(0)
> try(system("convert tmp/2fghp1353353270.ps tmp/2fghp1353353270.png",intern=TRUE))
character(0)
> try(system("convert tmp/3sjql1353353270.ps tmp/3sjql1353353270.png",intern=TRUE))
character(0)
> try(system("convert tmp/4wgqe1353353270.ps tmp/4wgqe1353353270.png",intern=TRUE))
character(0)
> try(system("convert tmp/5tiui1353353270.ps tmp/5tiui1353353270.png",intern=TRUE))
character(0)
> try(system("convert tmp/6b2zo1353353270.ps tmp/6b2zo1353353270.png",intern=TRUE))
character(0)
> try(system("convert tmp/7qeol1353353270.ps tmp/7qeol1353353270.png",intern=TRUE))
character(0)
> try(system("convert tmp/8527g1353353270.ps tmp/8527g1353353270.png",intern=TRUE))
character(0)
> try(system("convert tmp/9tme81353353270.ps tmp/9tme81353353270.png",intern=TRUE))
character(0)
> try(system("convert tmp/10bc9e1353353270.ps tmp/10bc9e1353353270.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
5.851 1.107 6.974