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
+ ,78051
+ ,47645
+ ,15545
+ ,35668
+ ,575093
+ ,129842
+ ,301647
+ ,5234
+ ,75481
+ ,45970
+ ,15001
+ ,35589
+ ,557560
+ ,129694
+ ,305353
+ ,5279
+ ,78926
+ ,48069
+ ,14961
+ ,35544
+ ,564478
+ ,130080
+ ,313665
+ ,5391
+ ,86241
+ ,53080
+ ,15245
+ ,35292
+ ,580523
+ ,131496
+ ,322402
+ ,5280
+ ,91993
+ ,57896
+ ,15656
+ ,35047
+ ,596594
+ ,131556
+ ,318280
+ ,5173
+ ,86452
+ ,54344
+ ,15577
+ ,34705
+ ,586570
+ ,128925
+ ,292852
+ ,4724
+ ,65271
+ ,40482
+ ,14630
+ ,34536
+ ,536214
+ ,127836
+ ,287481
+ ,4554
+ ,60348
+ ,37110
+ ,14336
+ ,33596
+ ,523597
+ ,129164
+ ,295210
+ ,4713
+ ,63178
+ ,39263
+ ,14834
+ ,34149
+ ,536535
+ ,129531
+ ,295650
+ ,4811
+ ,62653
+ ,38889
+ ,14921
+ ,33567
+ ,536322
+ ,128548
+ ,292919
+ ,4668
+ ,63583
+ ,39593
+ ,14707
+ ,32881
+ ,532638
+ ,127330
+ ,290649
+ ,4516
+ ,63376
+ ,39305
+ ,14516
+ ,32351
+ ,528222
+ ,123815
+ ,281687
+ ,4203
+ ,65008
+ ,40560
+ ,14055
+ ,31576
+ ,516141
+ ,124393
+ ,270336
+ ,4016
+ ,62100
+ ,38306
+ ,13493
+ ,31544
+ ,501866
+ ,123707
+ ,271420
+ ,3993
+ ,65936
+ ,40911
+ ,13528
+ ,31583
+ ,506174
+ ,123736
+ ,278183
+ ,3971
+ ,71621
+ ,44700
+ ,13719
+ ,30686
+ ,517945
+ ,124507
+ ,284913
+ ,3838
+ ,78903
+ ,50328
+ ,14170
+ ,31097
+ ,533590
+ ,125005
+ ,283487
+ ,3891
+ ,74755
+ ,47499
+ ,14009
+ ,31123
+ ,528379
+ ,121383
+ ,256677
+ ,3306
+ ,55511
+ ,34446
+ ,13159
+ ,30850
+ ,477580
+ ,121200
+ ,252945
+ ,3235
+ ,51888
+ ,31434
+ ,12927
+ ,30397
+ ,469357
+ ,125249
+ ,264963
+ ,3404
+ ,55738
+ ,34066
+ ,13510
+ ,30783
+ ,490243
+ ,125253
+ ,265988
+ ,3400
+ ,57261
+ ,35044
+ ,13520
+ ,30600
+ ,492622
+ ,127977
+ ,274857
+ ,3447
+ ,60086
+ ,37040
+ ,14089
+ ,30552
+ ,507561
+ ,128984
+ ,279650
+ ,3431
+ ,63070
+ ,38706
+ ,14251
+ ,30967
+ ,516922
+ ,126770
+ ,276715
+ ,3321
+ ,66061
+ ,40430
+ ,13980
+ ,30732
+ ,514258
+ ,126448
+ ,273887
+ ,3189
+ ,64973
+ ,39613
+ ,13715
+ ,30823
+ ,509846
+ ,127845
+ ,282308
+ ,3256
+ ,71770
+ ,44236
+ ,14112
+ ,31035
+ ,527070
+ ,128818
+ ,289847
+ ,3290
+ ,77712
+ ,47859
+ ,14289
+ ,30991
+ ,541657
+ ,132127
+ ,301101
+ ,3475
+ ,85265
+ ,53711
+ ,15020
+ ,31078
+ ,564591
+ ,132338
+ ,297008
+ ,3454
+ ,80140
+ ,50352
+ ,14860
+ ,31016
+ ,555362
+ ,126645
+ ,268909
+ ,2806
+ ,58921
+ ,36142
+ ,13800
+ ,30387
+ ,498662
+ ,130625
+ ,278383
+ ,2777
+ ,57395
+ ,34819
+ ,14431
+ ,30204
+ ,511038
+ ,133506
+ ,286226
+ ,2865
+ ,60925
+ ,37353
+ ,14944
+ ,30318
+ ,525919
+ ,135277
+ ,288936
+ ,2924
+ ,61682
+ ,37550
+ ,15083
+ ,30695
+ ,531673
+ ,137664
+ ,298953
+ ,3011
+ ,66161
+ ,40462
+ ,15707
+ ,30369
+ ,548854
+ ,139821
+ ,305837
+ ,3099
+ ,68713
+ ,41753
+ ,15954
+ ,30251
+ ,560576
+ ,138440
+ ,301979
+ ,2988
+ ,71442
+ ,43437
+ ,15631
+ ,29782
+ ,557274
+ ,139879
+ ,306281
+ ,3032
+ ,73898
+ ,44784
+ ,15813
+ ,29871
+ ,565742
+ ,142256
+ ,317057
+ ,3131
+ ,81482
+ ,49537
+ ,16356
+ ,30474
+ ,587625
+ ,146322
+ ,334780
+ ,3343
+ ,90533
+ ,54974
+ ,17086
+ ,31195
+ ,619916
+ ,146389
+ ,335895
+ ,3275
+ ,94794
+ ,58535
+ ,17302
+ ,31429
+ ,625809
+ ,147841
+ ,333874
+ ,3243
+ ,88780
+ ,54762
+ ,17247
+ ,31825
+ ,619567
+ ,146449
+ ,311028
+ ,2897
+ ,67281
+ ,40738
+ ,16398
+ ,31786
+ ,572942
+ ,147960
+ ,311767
+ ,2818
+ ,63724
+ ,38052
+ ,16590
+ ,32734
+ ,572775
+ ,148487
+ ,312575
+ ,2836
+ ,64361
+ ,38436
+ ,16673
+ ,32109
+ ,574205
+ ,149802
+ ,315040
+ ,2721
+ ,65465
+ ,36993
+ ,16962
+ ,32530
+ ,579799
+ ,151387
+ ,320325
+ ,2742
+ ,68725
+ ,39056
+ ,17278
+ ,32357
+ ,590072
+ ,151936
+ ,321178
+ ,2707
+ ,70782
+ ,39996
+ ,17224
+ ,32288
+ ,593408)
+ ,dim=c(8
+ ,48)
+ ,dimnames=list(c('Basisonderwijs(lager_1ste_graad_secundair)'
+ ,'Secundair_onderwijs(2de
+ ,3de
+ ,4de_graad)'
+ ,'Duaal_onderwijs'
+ ,'Hoger_onderwijs'
+ ,'Hoger_onderwijs(Bachelor)'
+ ,'Leercontract'
+ ,'Andere_studies'
+ ,'Werkloosheid_totaal')
+ ,1:48))
> y <- array(NA,dim=c(8,48),dimnames=list(c('Basisonderwijs(lager_1ste_graad_secundair)','Secundair_onderwijs(2de,3de,4de_graad)','Duaal_onderwijs','Hoger_onderwijs','Hoger_onderwijs(Bachelor)','Leercontract','Andere_studies','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 = '8'
> 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 Hoger_onderwijs
1 312991 5599 78051
2 301647 5234 75481
3 305353 5279 78926
4 313665 5391 86241
5 322402 5280 91993
6 318280 5173 86452
7 292852 4724 65271
8 287481 4554 60348
9 295210 4713 63178
10 295650 4811 62653
11 292919 4668 63583
12 290649 4516 63376
13 281687 4203 65008
14 270336 4016 62100
15 271420 3993 65936
16 278183 3971 71621
17 284913 3838 78903
18 283487 3891 74755
19 256677 3306 55511
20 252945 3235 51888
21 264963 3404 55738
22 265988 3400 57261
23 274857 3447 60086
24 279650 3431 63070
25 276715 3321 66061
26 273887 3189 64973
27 282308 3256 71770
28 289847 3290 77712
29 301101 3475 85265
30 297008 3454 80140
31 268909 2806 58921
32 278383 2777 57395
33 286226 2865 60925
34 288936 2924 61682
35 298953 3011 66161
36 305837 3099 68713
37 301979 2988 71442
38 306281 3032 73898
39 317057 3131 81482
40 334780 3343 90533
41 335895 3275 94794
42 333874 3243 88780
43 311028 2897 67281
44 311767 2818 63724
45 312575 2836 64361
46 315040 2721 65465
47 320325 2742 68725
48 321178 2707 70782
Hoger_onderwijs(Bachelor) Leercontract Andere_studies
1 47645 15545 35668
2 45970 15001 35589
3 48069 14961 35544
4 53080 15245 35292
5 57896 15656 35047
6 54344 15577 34705
7 40482 14630 34536
8 37110 14336 33596
9 39263 14834 34149
10 38889 14921 33567
11 39593 14707 32881
12 39305 14516 32351
13 40560 14055 31576
14 38306 13493 31544
15 40911 13528 31583
16 44700 13719 30686
17 50328 14170 31097
18 47499 14009 31123
19 34446 13159 30850
20 31434 12927 30397
21 34066 13510 30783
22 35044 13520 30600
23 37040 14089 30552
24 38706 14251 30967
25 40430 13980 30732
26 39613 13715 30823
27 44236 14112 31035
28 47859 14289 30991
29 53711 15020 31078
30 50352 14860 31016
31 36142 13800 30387
32 34819 14431 30204
33 37353 14944 30318
34 37550 15083 30695
35 40462 15707 30369
36 41753 15954 30251
37 43437 15631 29782
38 44784 15813 29871
39 49537 16356 30474
40 54974 17086 31195
41 58535 17302 31429
42 54762 17247 31825
43 40738 16398 31786
44 38052 16590 32734
45 38436 16673 32109
46 36993 16962 32530
47 39056 17278 32357
48 39996 17224 32288
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept)
3.079e-11
`Basisonderwijs(lager_1ste_graad_secundair)`
1.000e+00
`Secundair_onderwijs(2de,3de,4de_graad)`
1.000e+00
Duaal_onderwijs
4.377e-14
Hoger_onderwijs
1.000e+00
`Hoger_onderwijs(Bachelor)`
-2.545e-15
Leercontract
1.000e+00
Andere_studies
1.000e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3.098e-11 -6.729e-12 -1.832e-12 2.484e-12 1.171e-10
Coefficients:
Estimate Std. Error t value
(Intercept) 3.079e-11 1.826e-10 1.690e-01
`Basisonderwijs(lager_1ste_graad_secundair)` 1.000e+00 3.856e-15 2.594e+14
`Secundair_onderwijs(2de,3de,4de_graad)` 1.000e+00 2.218e-15 4.508e+14
Duaal_onderwijs 4.377e-14 2.045e-14 2.140e+00
Hoger_onderwijs 1.000e+00 3.714e-15 2.692e+14
`Hoger_onderwijs(Bachelor)` -2.545e-15 5.885e-15 -4.320e-01
Leercontract 1.000e+00 3.831e-14 2.610e+13
Andere_studies 1.000e+00 5.823e-15 1.717e+14
Pr(>|t|)
(Intercept) 0.8670
`Basisonderwijs(lager_1ste_graad_secundair)` <2e-16 ***
`Secundair_onderwijs(2de,3de,4de_graad)` <2e-16 ***
Duaal_onderwijs 0.0385 *
Hoger_onderwijs <2e-16 ***
`Hoger_onderwijs(Bachelor)` 0.6677
Leercontract <2e-16 ***
Andere_studies <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.121e-11 on 40 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 2.126e+31 on 7 and 40 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,] 4.351772e-02 8.703545e-02 9.564823e-01
[2,] 8.576330e-01 2.847340e-01 1.423670e-01
[3,] 9.999999e-01 2.700845e-07 1.350422e-07
[4,] 6.555800e-01 6.888401e-01 3.444200e-01
[5,] 9.999847e-01 3.051884e-05 1.525942e-05
[6,] 9.999999e-01 2.007787e-07 1.003893e-07
[7,] 7.522972e-01 4.954057e-01 2.477028e-01
[8,] 8.965355e-02 1.793071e-01 9.103465e-01
[9,] 9.999990e-01 1.934589e-06 9.672946e-07
[10,] 1.000000e+00 1.909502e-08 9.547511e-09
[11,] 5.495507e-03 1.099101e-02 9.945045e-01
[12,] 3.236459e-04 6.472919e-04 9.996764e-01
[13,] 5.411440e-13 1.082288e-12 1.000000e+00
[14,] 5.899516e-05 1.179903e-04 9.999410e-01
[15,] 7.870109e-01 4.259782e-01 2.129891e-01
[16,] 4.391536e-03 8.783071e-03 9.956085e-01
[17,] 4.604331e-01 9.208663e-01 5.395669e-01
[18,] 4.351359e-01 8.702717e-01 5.648641e-01
[19,] 1.134552e-02 2.269103e-02 9.886545e-01
[20,] 1.949243e-01 3.898486e-01 8.050757e-01
[21,] 1.217661e-03 2.435322e-03 9.987823e-01
[22,] 9.838070e-01 3.238606e-02 1.619303e-02
[23,] 9.999624e-01 7.512329e-05 3.756165e-05
[24,] 7.935708e-01 4.128585e-01 2.064292e-01
[25,] 8.974972e-01 2.050055e-01 1.025028e-01
[26,] 7.316963e-01 5.366075e-01 2.683037e-01
[27,] 5.967696e-01 8.064608e-01 4.032304e-01
> postscript(file="/var/fisher/rcomp/tmp/18pci1353349734.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/fisher/rcomp/tmp/2ovka1353349734.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/fisher/rcomp/tmp/3akaj1353349734.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/fisher/rcomp/tmp/4kkmr1353349734.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/fisher/rcomp/tmp/5be8i1353349734.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
1.170547e-10 -3.098400e-11 -1.428165e-11 -1.630650e-11 -1.437349e-11
6 7 8 9 10
-1.513666e-11 -6.259773e-12 -3.160587e-12 -3.359262e-12 -1.336082e-11
11 12 13 14 15
-1.046276e-11 -4.052549e-12 5.602780e-14 -1.577355e-11 -1.243953e-11
16 17 18 19 20
-3.952677e-12 2.190257e-12 1.683306e-12 -1.818244e-12 -5.742636e-12
21 22 23 24 25
-1.618374e-12 -4.882997e-12 -8.138566e-12 9.699249e-13 2.627181e-12
26 27 28 29 30
8.170460e-12 1.356480e-11 1.469271e-11 5.288731e-12 2.436109e-12
31 32 33 34 35
1.382232e-11 1.458919e-11 1.252284e-11 8.541233e-12 2.637130e-12
36 37 38 39 40
-1.453627e-12 -1.845403e-12 -2.367700e-12 -3.467314e-12 -3.383714e-12
41 42 43 44 45
-1.630244e-12 -2.467825e-12 -8.773768e-13 6.797475e-12 -1.119087e-12
46 47 48
1.757458e-13 -1.144057e-11 -1.166266e-11
> postscript(file="/var/fisher/rcomp/tmp/6zsdq1353349734.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 1.170547e-10 NA
1 -3.098400e-11 1.170547e-10
2 -1.428165e-11 -3.098400e-11
3 -1.630650e-11 -1.428165e-11
4 -1.437349e-11 -1.630650e-11
5 -1.513666e-11 -1.437349e-11
6 -6.259773e-12 -1.513666e-11
7 -3.160587e-12 -6.259773e-12
8 -3.359262e-12 -3.160587e-12
9 -1.336082e-11 -3.359262e-12
10 -1.046276e-11 -1.336082e-11
11 -4.052549e-12 -1.046276e-11
12 5.602780e-14 -4.052549e-12
13 -1.577355e-11 5.602780e-14
14 -1.243953e-11 -1.577355e-11
15 -3.952677e-12 -1.243953e-11
16 2.190257e-12 -3.952677e-12
17 1.683306e-12 2.190257e-12
18 -1.818244e-12 1.683306e-12
19 -5.742636e-12 -1.818244e-12
20 -1.618374e-12 -5.742636e-12
21 -4.882997e-12 -1.618374e-12
22 -8.138566e-12 -4.882997e-12
23 9.699249e-13 -8.138566e-12
24 2.627181e-12 9.699249e-13
25 8.170460e-12 2.627181e-12
26 1.356480e-11 8.170460e-12
27 1.469271e-11 1.356480e-11
28 5.288731e-12 1.469271e-11
29 2.436109e-12 5.288731e-12
30 1.382232e-11 2.436109e-12
31 1.458919e-11 1.382232e-11
32 1.252284e-11 1.458919e-11
33 8.541233e-12 1.252284e-11
34 2.637130e-12 8.541233e-12
35 -1.453627e-12 2.637130e-12
36 -1.845403e-12 -1.453627e-12
37 -2.367700e-12 -1.845403e-12
38 -3.467314e-12 -2.367700e-12
39 -3.383714e-12 -3.467314e-12
40 -1.630244e-12 -3.383714e-12
41 -2.467825e-12 -1.630244e-12
42 -8.773768e-13 -2.467825e-12
43 6.797475e-12 -8.773768e-13
44 -1.119087e-12 6.797475e-12
45 1.757458e-13 -1.119087e-12
46 -1.144057e-11 1.757458e-13
47 -1.166266e-11 -1.144057e-11
48 NA -1.166266e-11
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -3.098400e-11 1.170547e-10
[2,] -1.428165e-11 -3.098400e-11
[3,] -1.630650e-11 -1.428165e-11
[4,] -1.437349e-11 -1.630650e-11
[5,] -1.513666e-11 -1.437349e-11
[6,] -6.259773e-12 -1.513666e-11
[7,] -3.160587e-12 -6.259773e-12
[8,] -3.359262e-12 -3.160587e-12
[9,] -1.336082e-11 -3.359262e-12
[10,] -1.046276e-11 -1.336082e-11
[11,] -4.052549e-12 -1.046276e-11
[12,] 5.602780e-14 -4.052549e-12
[13,] -1.577355e-11 5.602780e-14
[14,] -1.243953e-11 -1.577355e-11
[15,] -3.952677e-12 -1.243953e-11
[16,] 2.190257e-12 -3.952677e-12
[17,] 1.683306e-12 2.190257e-12
[18,] -1.818244e-12 1.683306e-12
[19,] -5.742636e-12 -1.818244e-12
[20,] -1.618374e-12 -5.742636e-12
[21,] -4.882997e-12 -1.618374e-12
[22,] -8.138566e-12 -4.882997e-12
[23,] 9.699249e-13 -8.138566e-12
[24,] 2.627181e-12 9.699249e-13
[25,] 8.170460e-12 2.627181e-12
[26,] 1.356480e-11 8.170460e-12
[27,] 1.469271e-11 1.356480e-11
[28,] 5.288731e-12 1.469271e-11
[29,] 2.436109e-12 5.288731e-12
[30,] 1.382232e-11 2.436109e-12
[31,] 1.458919e-11 1.382232e-11
[32,] 1.252284e-11 1.458919e-11
[33,] 8.541233e-12 1.252284e-11
[34,] 2.637130e-12 8.541233e-12
[35,] -1.453627e-12 2.637130e-12
[36,] -1.845403e-12 -1.453627e-12
[37,] -2.367700e-12 -1.845403e-12
[38,] -3.467314e-12 -2.367700e-12
[39,] -3.383714e-12 -3.467314e-12
[40,] -1.630244e-12 -3.383714e-12
[41,] -2.467825e-12 -1.630244e-12
[42,] -8.773768e-13 -2.467825e-12
[43,] 6.797475e-12 -8.773768e-13
[44,] -1.119087e-12 6.797475e-12
[45,] 1.757458e-13 -1.119087e-12
[46,] -1.144057e-11 1.757458e-13
[47,] -1.166266e-11 -1.144057e-11
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -3.098400e-11 1.170547e-10
2 -1.428165e-11 -3.098400e-11
3 -1.630650e-11 -1.428165e-11
4 -1.437349e-11 -1.630650e-11
5 -1.513666e-11 -1.437349e-11
6 -6.259773e-12 -1.513666e-11
7 -3.160587e-12 -6.259773e-12
8 -3.359262e-12 -3.160587e-12
9 -1.336082e-11 -3.359262e-12
10 -1.046276e-11 -1.336082e-11
11 -4.052549e-12 -1.046276e-11
12 5.602780e-14 -4.052549e-12
13 -1.577355e-11 5.602780e-14
14 -1.243953e-11 -1.577355e-11
15 -3.952677e-12 -1.243953e-11
16 2.190257e-12 -3.952677e-12
17 1.683306e-12 2.190257e-12
18 -1.818244e-12 1.683306e-12
19 -5.742636e-12 -1.818244e-12
20 -1.618374e-12 -5.742636e-12
21 -4.882997e-12 -1.618374e-12
22 -8.138566e-12 -4.882997e-12
23 9.699249e-13 -8.138566e-12
24 2.627181e-12 9.699249e-13
25 8.170460e-12 2.627181e-12
26 1.356480e-11 8.170460e-12
27 1.469271e-11 1.356480e-11
28 5.288731e-12 1.469271e-11
29 2.436109e-12 5.288731e-12
30 1.382232e-11 2.436109e-12
31 1.458919e-11 1.382232e-11
32 1.252284e-11 1.458919e-11
33 8.541233e-12 1.252284e-11
34 2.637130e-12 8.541233e-12
35 -1.453627e-12 2.637130e-12
36 -1.845403e-12 -1.453627e-12
37 -2.367700e-12 -1.845403e-12
38 -3.467314e-12 -2.367700e-12
39 -3.383714e-12 -3.467314e-12
40 -1.630244e-12 -3.383714e-12
41 -2.467825e-12 -1.630244e-12
42 -8.773768e-13 -2.467825e-12
43 6.797475e-12 -8.773768e-13
44 -1.119087e-12 6.797475e-12
45 1.757458e-13 -1.119087e-12
46 -1.144057e-11 1.757458e-13
47 -1.166266e-11 -1.144057e-11
> 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/fisher/rcomp/tmp/75lor1353349734.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/fisher/rcomp/tmp/8hh7l1353349734.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/fisher/rcomp/tmp/9k6o61353349734.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/fisher/rcomp/tmp/10b8881353349734.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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/fisher/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/fisher/rcomp/tmp/116zh61353349734.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/fisher/rcomp/tmp/12il031353349734.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/fisher/rcomp/tmp/13jkr51353349734.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/fisher/rcomp/tmp/14z8c21353349734.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/fisher/rcomp/tmp/157y051353349734.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/fisher/rcomp/tmp/16c7kw1353349734.tab")
+ }
>
> try(system("convert tmp/18pci1353349734.ps tmp/18pci1353349734.png",intern=TRUE))
character(0)
> try(system("convert tmp/2ovka1353349734.ps tmp/2ovka1353349734.png",intern=TRUE))
character(0)
> try(system("convert tmp/3akaj1353349734.ps tmp/3akaj1353349734.png",intern=TRUE))
character(0)
> try(system("convert tmp/4kkmr1353349734.ps tmp/4kkmr1353349734.png",intern=TRUE))
character(0)
> try(system("convert tmp/5be8i1353349734.ps tmp/5be8i1353349734.png",intern=TRUE))
character(0)
> try(system("convert tmp/6zsdq1353349734.ps tmp/6zsdq1353349734.png",intern=TRUE))
character(0)
> try(system("convert tmp/75lor1353349734.ps tmp/75lor1353349734.png",intern=TRUE))
character(0)
> try(system("convert tmp/8hh7l1353349734.ps tmp/8hh7l1353349734.png",intern=TRUE))
character(0)
> try(system("convert tmp/9k6o61353349734.ps tmp/9k6o61353349734.png",intern=TRUE))
character(0)
> try(system("convert tmp/10b8881353349734.ps tmp/10b8881353349734.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
5.689 1.294 6.979