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
+ ,35668
+ ,575093
+ ,129842
+ ,301647
+ ,5234
+ ,45970
+ ,15001
+ ,35589
+ ,557560
+ ,129694
+ ,305353
+ ,5279
+ ,48069
+ ,14961
+ ,35544
+ ,564478
+ ,130080
+ ,313665
+ ,5391
+ ,53080
+ ,15245
+ ,35292
+ ,580523
+ ,131496
+ ,322402
+ ,5280
+ ,57896
+ ,15656
+ ,35047
+ ,596594
+ ,131556
+ ,318280
+ ,5173
+ ,54344
+ ,15577
+ ,34705
+ ,586570
+ ,128925
+ ,292852
+ ,4724
+ ,40482
+ ,14630
+ ,34536
+ ,536214
+ ,127836
+ ,287481
+ ,4554
+ ,37110
+ ,14336
+ ,33596
+ ,523597
+ ,129164
+ ,295210
+ ,4713
+ ,39263
+ ,14834
+ ,34149
+ ,536535
+ ,129531
+ ,295650
+ ,4811
+ ,38889
+ ,14921
+ ,33567
+ ,536322
+ ,128548
+ ,292919
+ ,4668
+ ,39593
+ ,14707
+ ,32881
+ ,532638
+ ,127330
+ ,290649
+ ,4516
+ ,39305
+ ,14516
+ ,32351
+ ,528222
+ ,123815
+ ,281687
+ ,4203
+ ,40560
+ ,14055
+ ,31576
+ ,516141
+ ,124393
+ ,270336
+ ,4016
+ ,38306
+ ,13493
+ ,31544
+ ,501866
+ ,123707
+ ,271420
+ ,3993
+ ,40911
+ ,13528
+ ,31583
+ ,506174
+ ,123736
+ ,278183
+ ,3971
+ ,44700
+ ,13719
+ ,30686
+ ,517945
+ ,124507
+ ,284913
+ ,3838
+ ,50328
+ ,14170
+ ,31097
+ ,533590
+ ,125005
+ ,283487
+ ,3891
+ ,47499
+ ,14009
+ ,31123
+ ,528379
+ ,121383
+ ,256677
+ ,3306
+ ,34446
+ ,13159
+ ,30850
+ ,477580
+ ,121200
+ ,252945
+ ,3235
+ ,31434
+ ,12927
+ ,30397
+ ,469357
+ ,125249
+ ,264963
+ ,3404
+ ,34066
+ ,13510
+ ,30783
+ ,490243
+ ,125253
+ ,265988
+ ,3400
+ ,35044
+ ,13520
+ ,30600
+ ,492622
+ ,127977
+ ,274857
+ ,3447
+ ,37040
+ ,14089
+ ,30552
+ ,507561
+ ,128984
+ ,279650
+ ,3431
+ ,38706
+ ,14251
+ ,30967
+ ,516922
+ ,126770
+ ,276715
+ ,3321
+ ,40430
+ ,13980
+ ,30732
+ ,514258
+ ,126448
+ ,273887
+ ,3189
+ ,39613
+ ,13715
+ ,30823
+ ,509846
+ ,127845
+ ,282308
+ ,3256
+ ,44236
+ ,14112
+ ,31035
+ ,527070
+ ,128818
+ ,289847
+ ,3290
+ ,47859
+ ,14289
+ ,30991
+ ,541657
+ ,132127
+ ,301101
+ ,3475
+ ,53711
+ ,15020
+ ,31078
+ ,564591
+ ,132338
+ ,297008
+ ,3454
+ ,50352
+ ,14860
+ ,31016
+ ,555362
+ ,126645
+ ,268909
+ ,2806
+ ,36142
+ ,13800
+ ,30387
+ ,498662
+ ,130625
+ ,278383
+ ,2777
+ ,34819
+ ,14431
+ ,30204
+ ,511038
+ ,133506
+ ,286226
+ ,2865
+ ,37353
+ ,14944
+ ,30318
+ ,525919
+ ,135277
+ ,288936
+ ,2924
+ ,37550
+ ,15083
+ ,30695
+ ,531673
+ ,137664
+ ,298953
+ ,3011
+ ,40462
+ ,15707
+ ,30369
+ ,548854
+ ,139821
+ ,305837
+ ,3099
+ ,41753
+ ,15954
+ ,30251
+ ,560576
+ ,138440
+ ,301979
+ ,2988
+ ,43437
+ ,15631
+ ,29782
+ ,557274
+ ,139879
+ ,306281
+ ,3032
+ ,44784
+ ,15813
+ ,29871
+ ,565742
+ ,142256
+ ,317057
+ ,3131
+ ,49537
+ ,16356
+ ,30474
+ ,587625
+ ,146322
+ ,334780
+ ,3343
+ ,54974
+ ,17086
+ ,31195
+ ,619916
+ ,146389
+ ,335895
+ ,3275
+ ,58535
+ ,17302
+ ,31429
+ ,625809
+ ,147841
+ ,333874
+ ,3243
+ ,54762
+ ,17247
+ ,31825
+ ,619567
+ ,146449
+ ,311028
+ ,2897
+ ,40738
+ ,16398
+ ,31786
+ ,572942
+ ,147960
+ ,311767
+ ,2818
+ ,38052
+ ,16590
+ ,32734
+ ,572775
+ ,148487
+ ,312575
+ ,2836
+ ,38436
+ ,16673
+ ,32109
+ ,574205
+ ,149802
+ ,315040
+ ,2721
+ ,36993
+ ,16962
+ ,32530
+ ,579799
+ ,151387
+ ,320325
+ ,2742
+ ,39056
+ ,17278
+ ,32357
+ ,590072
+ ,151936
+ ,321178
+ ,2707
+ ,39996
+ ,17224
+ ,32288
+ ,593408)
+ ,dim=c(7
+ ,48)
+ ,dimnames=list(c('Basisonderwijs(lager_1ste_graad_secundair)'
+ ,'Secundair_onderwijs(2de
+ ,3de
+ ,4de_graad)'
+ ,'Duaal_onderwijs'
+ ,'Hoger_onderwijs(Bachelor)'
+ ,'Leercontract'
+ ,'Andere_studies'
+ ,'Werkloosheid_totaal')
+ ,1:48))
> y <- array(NA,dim=c(7,48),dimnames=list(c('Basisonderwijs(lager_1ste_graad_secundair)','Secundair_onderwijs(2de,3de,4de_graad)','Duaal_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 = '7'
> 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 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)
-2.808e+04
`Basisonderwijs(lager_1ste_graad_secundair)`
1.309e+00
`Secundair_onderwijs(2de,3de,4de_graad)`
1.023e+00
Duaal_onderwijs
-2.940e-01
`Hoger_onderwijs(Bachelor)`
1.510e+00
Leercontract
-6.437e-01
Andere_studies
1.348e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1907.34 -521.02 8.78 580.64 1932.63
Coefficients:
Estimate Std. Error t value
(Intercept) -2.808e+04 6.304e+03 -4.453
`Basisonderwijs(lager_1ste_graad_secundair)` 1.309e+00 1.547e-01 8.462
`Secundair_onderwijs(2de,3de,4de_graad)` 1.023e+00 9.320e-02 10.977
Duaal_onderwijs -2.940e-01 8.587e-01 -0.342
`Hoger_onderwijs(Bachelor)` 1.510e+00 7.478e-02 20.196
Leercontract -6.437e-01 1.590e+00 -0.405
Andere_studies 1.348e+00 2.387e-01 5.646
Pr(>|t|)
(Intercept) 6.38e-05 ***
`Basisonderwijs(lager_1ste_graad_secundair)` 1.55e-10 ***
`Secundair_onderwijs(2de,3de,4de_graad)` 8.88e-14 ***
Duaal_onderwijs 0.734
`Hoger_onderwijs(Bachelor)` < 2e-16 ***
Leercontract 0.688
Andere_studies 1.38e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 891.8 on 41 degrees of freedom
Multiple R-squared: 0.9995, Adjusted R-squared: 0.9994
F-statistic: 1.402e+04 on 6 and 41 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,] 1.375464e-02 2.750927e-02 0.98624536
[2,] 2.431559e-03 4.863119e-03 0.99756844
[3,] 1.325874e-03 2.651749e-03 0.99867413
[4,] 2.696767e-04 5.393534e-04 0.99973032
[5,] 6.153645e-05 1.230729e-04 0.99993846
[6,] 1.024845e-05 2.049689e-05 0.99998975
[7,] 1.152789e-05 2.305579e-05 0.99998847
[8,] 1.319336e-05 2.638672e-05 0.99998681
[9,] 3.026645e-06 6.053290e-06 0.99999697
[10,] 2.490113e-05 4.980225e-05 0.99997510
[11,] 1.888055e-03 3.776109e-03 0.99811195
[12,] 1.074409e-02 2.148817e-02 0.98925591
[13,] 1.089913e-02 2.179827e-02 0.98910087
[14,] 9.113003e-03 1.822601e-02 0.99088700
[15,] 7.838925e-03 1.567785e-02 0.99216107
[16,] 1.725347e-02 3.450694e-02 0.98274653
[17,] 1.525593e-02 3.051186e-02 0.98474407
[18,] 1.181865e-02 2.363731e-02 0.98818135
[19,] 9.799426e-03 1.959885e-02 0.99020057
[20,] 8.754561e-03 1.750912e-02 0.99124544
[21,] 1.974884e-01 3.949768e-01 0.80251158
[22,] 6.056264e-01 7.887472e-01 0.39437362
[23,] 6.664828e-01 6.670345e-01 0.33351725
[24,] 8.744462e-01 2.511076e-01 0.12555382
[25,] 9.692275e-01 6.154494e-02 0.03077247
[26,] 9.499393e-01 1.001213e-01 0.05006066
[27,] 9.029871e-01 1.940258e-01 0.09701290
[28,] 9.107731e-01 1.784537e-01 0.08922685
[29,] 9.130901e-01 1.738198e-01 0.08690989
> postscript(file="/var/fisher/rcomp/tmp/1wu6l1353352667.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/2l7h71353352667.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/3tsnw1353352667.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/4vr9l1353352667.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/5sixo1353352667.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 6
663.088435 836.585162 1035.356846 1059.459798 -372.706953 -515.546203
7 8 9 10 11 12
-992.865089 -569.208330 -906.368503 -615.990897 -537.438988 -54.703835
13 14 15 16 17 18
395.447379 5.929208 132.343252 550.096755 -501.525723 -756.499839
19 20 21 22 23 24
-24.540199 799.033084 18.589832 118.690730 -151.773275 11.626163
25 26 27 28 29 30
755.527530 559.992718 347.827604 660.099451 -681.779019 -952.723984
31 32 33 34 35 36
-17.791173 97.144736 -442.163903 -478.316193 -201.401712 48.067206
37 38 39 40 41 42
349.430668 508.091341 642.599063 828.316573 -81.331920 -1037.637786
43 44 45 46 47 48
-1884.908865 -1907.343985 -1672.838555 1441.629258 1559.802344 1932.629790
> postscript(file="/var/fisher/rcomp/tmp/6lrac1353352667.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 663.088435 NA
1 836.585162 663.088435
2 1035.356846 836.585162
3 1059.459798 1035.356846
4 -372.706953 1059.459798
5 -515.546203 -372.706953
6 -992.865089 -515.546203
7 -569.208330 -992.865089
8 -906.368503 -569.208330
9 -615.990897 -906.368503
10 -537.438988 -615.990897
11 -54.703835 -537.438988
12 395.447379 -54.703835
13 5.929208 395.447379
14 132.343252 5.929208
15 550.096755 132.343252
16 -501.525723 550.096755
17 -756.499839 -501.525723
18 -24.540199 -756.499839
19 799.033084 -24.540199
20 18.589832 799.033084
21 118.690730 18.589832
22 -151.773275 118.690730
23 11.626163 -151.773275
24 755.527530 11.626163
25 559.992718 755.527530
26 347.827604 559.992718
27 660.099451 347.827604
28 -681.779019 660.099451
29 -952.723984 -681.779019
30 -17.791173 -952.723984
31 97.144736 -17.791173
32 -442.163903 97.144736
33 -478.316193 -442.163903
34 -201.401712 -478.316193
35 48.067206 -201.401712
36 349.430668 48.067206
37 508.091341 349.430668
38 642.599063 508.091341
39 828.316573 642.599063
40 -81.331920 828.316573
41 -1037.637786 -81.331920
42 -1884.908865 -1037.637786
43 -1907.343985 -1884.908865
44 -1672.838555 -1907.343985
45 1441.629258 -1672.838555
46 1559.802344 1441.629258
47 1932.629790 1559.802344
48 NA 1932.629790
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 836.585162 663.088435
[2,] 1035.356846 836.585162
[3,] 1059.459798 1035.356846
[4,] -372.706953 1059.459798
[5,] -515.546203 -372.706953
[6,] -992.865089 -515.546203
[7,] -569.208330 -992.865089
[8,] -906.368503 -569.208330
[9,] -615.990897 -906.368503
[10,] -537.438988 -615.990897
[11,] -54.703835 -537.438988
[12,] 395.447379 -54.703835
[13,] 5.929208 395.447379
[14,] 132.343252 5.929208
[15,] 550.096755 132.343252
[16,] -501.525723 550.096755
[17,] -756.499839 -501.525723
[18,] -24.540199 -756.499839
[19,] 799.033084 -24.540199
[20,] 18.589832 799.033084
[21,] 118.690730 18.589832
[22,] -151.773275 118.690730
[23,] 11.626163 -151.773275
[24,] 755.527530 11.626163
[25,] 559.992718 755.527530
[26,] 347.827604 559.992718
[27,] 660.099451 347.827604
[28,] -681.779019 660.099451
[29,] -952.723984 -681.779019
[30,] -17.791173 -952.723984
[31,] 97.144736 -17.791173
[32,] -442.163903 97.144736
[33,] -478.316193 -442.163903
[34,] -201.401712 -478.316193
[35,] 48.067206 -201.401712
[36,] 349.430668 48.067206
[37,] 508.091341 349.430668
[38,] 642.599063 508.091341
[39,] 828.316573 642.599063
[40,] -81.331920 828.316573
[41,] -1037.637786 -81.331920
[42,] -1884.908865 -1037.637786
[43,] -1907.343985 -1884.908865
[44,] -1672.838555 -1907.343985
[45,] 1441.629258 -1672.838555
[46,] 1559.802344 1441.629258
[47,] 1932.629790 1559.802344
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 836.585162 663.088435
2 1035.356846 836.585162
3 1059.459798 1035.356846
4 -372.706953 1059.459798
5 -515.546203 -372.706953
6 -992.865089 -515.546203
7 -569.208330 -992.865089
8 -906.368503 -569.208330
9 -615.990897 -906.368503
10 -537.438988 -615.990897
11 -54.703835 -537.438988
12 395.447379 -54.703835
13 5.929208 395.447379
14 132.343252 5.929208
15 550.096755 132.343252
16 -501.525723 550.096755
17 -756.499839 -501.525723
18 -24.540199 -756.499839
19 799.033084 -24.540199
20 18.589832 799.033084
21 118.690730 18.589832
22 -151.773275 118.690730
23 11.626163 -151.773275
24 755.527530 11.626163
25 559.992718 755.527530
26 347.827604 559.992718
27 660.099451 347.827604
28 -681.779019 660.099451
29 -952.723984 -681.779019
30 -17.791173 -952.723984
31 97.144736 -17.791173
32 -442.163903 97.144736
33 -478.316193 -442.163903
34 -201.401712 -478.316193
35 48.067206 -201.401712
36 349.430668 48.067206
37 508.091341 349.430668
38 642.599063 508.091341
39 828.316573 642.599063
40 -81.331920 828.316573
41 -1037.637786 -81.331920
42 -1884.908865 -1037.637786
43 -1907.343985 -1884.908865
44 -1672.838555 -1907.343985
45 1441.629258 -1672.838555
46 1559.802344 1441.629258
47 1932.629790 1559.802344
> 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/70dq61353352667.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/8ps1p1353352667.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/9ef291353352667.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/107fp71353352667.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/11nz5c1353352667.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/12co0t1353352667.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/13g5qj1353352667.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/14i0ix1353352667.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/15f33o1353352667.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/165ynb1353352668.tab")
+ }
>
> try(system("convert tmp/1wu6l1353352667.ps tmp/1wu6l1353352667.png",intern=TRUE))
character(0)
> try(system("convert tmp/2l7h71353352667.ps tmp/2l7h71353352667.png",intern=TRUE))
character(0)
> try(system("convert tmp/3tsnw1353352667.ps tmp/3tsnw1353352667.png",intern=TRUE))
character(0)
> try(system("convert tmp/4vr9l1353352667.ps tmp/4vr9l1353352667.png",intern=TRUE))
character(0)
> try(system("convert tmp/5sixo1353352667.ps tmp/5sixo1353352667.png",intern=TRUE))
character(0)
> try(system("convert tmp/6lrac1353352667.ps tmp/6lrac1353352667.png",intern=TRUE))
character(0)
> try(system("convert tmp/70dq61353352667.ps tmp/70dq61353352667.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ps1p1353352667.ps tmp/8ps1p1353352667.png",intern=TRUE))
character(0)
> try(system("convert tmp/9ef291353352667.ps tmp/9ef291353352667.png",intern=TRUE))
character(0)
> try(system("convert tmp/107fp71353352667.ps tmp/107fp71353352667.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
5.726 1.296 7.036