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(01/2006
+ ,593408
+ ,151936
+ ,321178
+ ,489
+ ,39507
+ ,30786
+ ,02/2006
+ ,590072
+ ,151387
+ ,320325
+ ,495
+ ,38561
+ ,29669
+ ,03/2006
+ ,579799
+ ,149802
+ ,315040
+ ,494
+ ,36499
+ ,28472
+ ,04/2006
+ ,574205
+ ,148487
+ ,312575
+ ,550
+ ,37886
+ ,25925
+ ,05/2006
+ ,572775
+ ,147960
+ ,311767
+ ,612
+ ,37440
+ ,25672
+ ,06/2006
+ ,572942
+ ,146449
+ ,311028
+ ,662
+ ,40076
+ ,26543
+ ,07/2006
+ ,619567
+ ,147841
+ ,333874
+ ,808
+ ,53954
+ ,34018
+ ,08/2006
+ ,625809
+ ,146389
+ ,335895
+ ,885
+ ,57650
+ ,36259
+ ,09/2006
+ ,619916
+ ,146322
+ ,334780
+ ,973
+ ,54001
+ ,35559
+ ,10/2006
+ ,587625
+ ,142256
+ ,317057
+ ,974
+ ,48563
+ ,31945
+ ,11/2006
+ ,565742
+ ,139879
+ ,306281
+ ,949
+ ,43835
+ ,29114
+ ,12/2006
+ ,557274
+ ,138440
+ ,301979
+ ,949
+ ,42488
+ ,28005
+ ,01/2007
+ ,560576
+ ,139821
+ ,305837
+ ,951
+ ,40802
+ ,26960
+ ,02/2007
+ ,548854
+ ,137664
+ ,298953
+ ,986
+ ,39476
+ ,25699
+ ,03/2007
+ ,531673
+ ,135277
+ ,288936
+ ,945
+ ,36605
+ ,24132
+ ,04/2007
+ ,525919
+ ,133506
+ ,286226
+ ,945
+ ,36408
+ ,23572
+ ,05/2007
+ ,511038
+ ,130625
+ ,278383
+ ,917
+ ,33902
+ ,22576
+ ,06/2007
+ ,498662
+ ,126645
+ ,268909
+ ,982
+ ,35160
+ ,22779
+ ,07/2007
+ ,555362
+ ,132338
+ ,297008
+ ,1248
+ ,49104
+ ,29788
+ ,08/2007
+ ,564591
+ ,132127
+ ,301101
+ ,1438
+ ,52273
+ ,31554
+ ,09/2007
+ ,541657
+ ,128818
+ ,289847
+ ,1551
+ ,46308
+ ,29853
+ ,10/2007
+ ,527070
+ ,127845
+ ,282308
+ ,1517
+ ,42719
+ ,27534
+ ,11/2007
+ ,509846
+ ,126448
+ ,273887
+ ,1442
+ ,38171
+ ,25360
+ ,12/2007
+ ,514258
+ ,126770
+ ,276715
+ ,1418
+ ,39012
+ ,25631
+ ,01/2008
+ ,516922
+ ,128984
+ ,279650
+ ,1383
+ ,37323
+ ,24364
+ ,02/2008
+ ,507561
+ ,127977
+ ,274857
+ ,1354
+ ,35686
+ ,23046
+ ,03/2008
+ ,492622
+ ,125253
+ ,265988
+ ,1310
+ ,33734
+ ,22217
+ ,04/2008
+ ,490243
+ ,125249
+ ,264963
+ ,1269
+ ,32797
+ ,21672
+ ,05/2008
+ ,469357
+ ,121200
+ ,252945
+ ,1198
+ ,30236
+ ,20454
+ ,06/2008
+ ,477580
+ ,121383
+ ,256677
+ ,1257
+ ,33189
+ ,21065
+ ,07/2008
+ ,528379
+ ,125005
+ ,283487
+ ,1585
+ ,45914
+ ,27256
+ ,08/2008
+ ,533590
+ ,124507
+ ,284913
+ ,1662
+ ,48666
+ ,28575
+ ,09/2008
+ ,517945
+ ,123736
+ ,278183
+ ,1695
+ ,43005
+ ,26921
+ ,10/2008
+ ,506174
+ ,123707
+ ,271420
+ ,1610
+ ,39301
+ ,25025
+ ,11/2008
+ ,501866
+ ,124393
+ ,270336
+ ,1580
+ ,36726
+ ,23794
+ ,12/2008
+ ,516141
+ ,123815
+ ,281687
+ ,1584
+ ,38976
+ ,24448
+ ,01/2009
+ ,528222
+ ,127330
+ ,290649
+ ,1573
+ ,37732
+ ,24071
+ ,02/2009
+ ,532638
+ ,128548
+ ,292919
+ ,1633
+ ,37960
+ ,23990
+ ,03/2009
+ ,536322
+ ,129531
+ ,295650
+ ,1631
+ ,37258
+ ,23764
+ ,04/2009
+ ,536535
+ ,129164
+ ,295210
+ ,1652
+ ,37611
+ ,23915
+ ,05/2009
+ ,523597
+ ,127836
+ ,287481
+ ,1591
+ ,35519
+ ,23238
+ ,06/2009
+ ,536214
+ ,128925
+ ,292852
+ ,1652
+ ,38830
+ ,24789
+ ,07/2009
+ ,586570
+ ,131556
+ ,318280
+ ,2034
+ ,52310
+ ,32108
+ ,08/2009
+ ,596594
+ ,131496
+ ,322402
+ ,2266
+ ,55630
+ ,34097
+ ,09/2009
+ ,580523
+ ,130080
+ ,313665
+ ,2372
+ ,50708
+ ,33161
+ ,10/2009
+ ,564478
+ ,129694
+ ,305353
+ ,2237
+ ,45832
+ ,30857
+ ,11/2009
+ ,557560
+ ,129842
+ ,301647
+ ,2118
+ ,43852
+ ,29511
+ ,12/2009
+ ,575093
+ ,132838
+ ,312991
+ ,2150
+ ,45495
+ ,30406
+ ,01/2010
+ ,580112
+ ,147512
+ ,335839
+ ,2629
+ ,48300
+ ,29975
+ ,02/2010
+ ,574761
+ ,147292
+ ,332590
+ ,2584
+ ,47043
+ ,29504
+ ,03/2010
+ ,563250
+ ,146997
+ ,325896
+ ,2442
+ ,44032
+ ,28655
+ ,04/2010
+ ,551531
+ ,144952
+ ,318433
+ ,2383
+ ,42872
+ ,28129
+ ,05/2010
+ ,537034
+ ,142704
+ ,309351
+ ,2275
+ ,40866
+ ,27435
+ ,06/2010
+ ,544686
+ ,143288
+ ,312122
+ ,2368
+ ,43635
+ ,28881
+ ,07/2010
+ ,600901
+ ,147234
+ ,342116
+ ,2866
+ ,57022
+ ,36183
+ ,08/2010
+ ,604378
+ ,146713
+ ,342105
+ ,3084
+ ,59494
+ ,37516
+ ,09/2010
+ ,586111
+ ,144235
+ ,332239
+ ,3018
+ ,54715
+ ,37078
+ ,10/2010
+ ,563698
+ ,143059
+ ,320198
+ ,2805
+ ,49098
+ ,34251
+ ,11/2010
+ ,548604
+ ,141610
+ ,311980
+ ,2688
+ ,46251
+ ,32039
+ ,12/2010
+ ,551074
+ ,142279
+ ,313907
+ ,2658
+ ,45915
+ ,32081)
+ ,dim=c(7
+ ,60)
+ ,dimnames=list(c('month'
+ ,'Totaal_Belgie'
+ ,'Basisonderwijs'
+ ,'Secundair_onderwijs'
+ ,'Academische_bachelor'
+ ,'Professionele_bachelor'
+ ,'Master_doctoraat')
+ ,1:60))
> y <- array(NA,dim=c(7,60),dimnames=list(c('month','Totaal_Belgie','Basisonderwijs','Secundair_onderwijs','Academische_bachelor','Professionele_bachelor','Master_doctoraat'),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 = '2'
> 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
Totaal_Belgie month Basisonderwijs Secundair_onderwijs
1 593408 0.0004985045 151936 321178
2 590072 0.0009970090 151387 320325
3 579799 0.0014955135 149802 315040
4 574205 0.0019940179 148487 312575
5 572775 0.0024925224 147960 311767
6 572942 0.0029910269 146449 311028
7 619567 0.0034895314 147841 333874
8 625809 0.0039880359 146389 335895
9 619916 0.0044865404 146322 334780
10 587625 0.0049850449 142256 317057
11 565742 0.0054835494 139879 306281
12 557274 0.0059820538 138440 301979
13 560576 0.0004982561 139821 305837
14 548854 0.0009965122 137664 298953
15 531673 0.0014947683 135277 288936
16 525919 0.0019930244 133506 286226
17 511038 0.0024912805 130625 278383
18 498662 0.0029895366 126645 268909
19 555362 0.0034877927 132338 297008
20 564591 0.0039860488 132127 301101
21 541657 0.0044843049 128818 289847
22 527070 0.0049825610 127845 282308
23 509846 0.0054808171 126448 273887
24 514258 0.0059790732 126770 276715
25 516922 0.0004980080 128984 279650
26 507561 0.0009960159 127977 274857
27 492622 0.0014940239 125253 265988
28 490243 0.0019920319 125249 264963
29 469357 0.0024900398 121200 252945
30 477580 0.0029880478 121383 256677
31 528379 0.0034860558 125005 283487
32 533590 0.0039840637 124507 284913
33 517945 0.0044820717 123736 278183
34 506174 0.0049800797 123707 271420
35 501866 0.0054780876 124393 270336
36 516141 0.0059760956 123815 281687
37 528222 0.0004977601 127330 290649
38 532638 0.0009955202 128548 292919
39 536322 0.0014932802 129531 295650
40 536535 0.0019910403 129164 295210
41 523597 0.0024888004 127836 287481
42 536214 0.0029865605 128925 292852
43 586570 0.0034843206 131556 318280
44 596594 0.0039820806 131496 322402
45 580523 0.0044798407 130080 313665
46 564478 0.0049776008 129694 305353
47 557560 0.0054753609 129842 301647
48 575093 0.0059731210 132838 312991
49 580112 0.0004975124 147512 335839
50 574761 0.0009950249 147292 332590
51 563250 0.0014925373 146997 325896
52 551531 0.0019900498 144952 318433
53 537034 0.0024875622 142704 309351
54 544686 0.0029850746 143288 312122
55 600901 0.0034825871 147234 342116
56 604378 0.0039800995 146713 342105
57 586111 0.0044776119 144235 332239
58 563698 0.0049751244 143059 320198
59 548604 0.0054726368 141610 311980
60 551074 0.0059701493 142279 313907
Academische_bachelor Professionele_bachelor Master_doctoraat
1 489 39507 30786
2 495 38561 29669
3 494 36499 28472
4 550 37886 25925
5 612 37440 25672
6 662 40076 26543
7 808 53954 34018
8 885 57650 36259
9 973 54001 35559
10 974 48563 31945
11 949 43835 29114
12 949 42488 28005
13 951 40802 26960
14 986 39476 25699
15 945 36605 24132
16 945 36408 23572
17 917 33902 22576
18 982 35160 22779
19 1248 49104 29788
20 1438 52273 31554
21 1551 46308 29853
22 1517 42719 27534
23 1442 38171 25360
24 1418 39012 25631
25 1383 37323 24364
26 1354 35686 23046
27 1310 33734 22217
28 1269 32797 21672
29 1198 30236 20454
30 1257 33189 21065
31 1585 45914 27256
32 1662 48666 28575
33 1695 43005 26921
34 1610 39301 25025
35 1580 36726 23794
36 1584 38976 24448
37 1573 37732 24071
38 1633 37960 23990
39 1631 37258 23764
40 1652 37611 23915
41 1591 35519 23238
42 1652 38830 24789
43 2034 52310 32108
44 2266 55630 34097
45 2372 50708 33161
46 2237 45832 30857
47 2118 43852 29511
48 2150 45495 30406
49 2629 48300 29975
50 2584 47043 29504
51 2442 44032 28655
52 2383 42872 28129
53 2275 40866 27435
54 2368 43635 28881
55 2866 57022 36183
56 3084 59494 37516
57 3018 54715 37078
58 2805 49098 34251
59 2688 46251 32039
60 2658 45915 32081
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) month Basisonderwijs
1.506e+05 1.414e+05 -1.508e+00
Secundair_onderwijs Academische_bachelor Professionele_bachelor
1.914e+00 -1.935e+01 -4.849e-01
Master_doctoraat
2.782e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-6245 -1653 -76 1497 6749
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.506e+05 8.120e+03 18.547 < 2e-16 ***
month 1.414e+05 3.061e+05 0.462 0.6459
Basisonderwijs -1.508e+00 1.294e-01 -11.650 3.09e-16 ***
Secundair_onderwijs 1.914e+00 7.037e-02 27.201 < 2e-16 ***
Academische_bachelor -1.935e+01 6.830e-01 -28.335 < 2e-16 ***
Professionele_bachelor -4.849e-01 2.113e-01 -2.295 0.0257 *
Master_doctoraat 2.782e+00 3.944e-01 7.053 3.71e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2884 on 53 degrees of freedom
Multiple R-squared: 0.9943, Adjusted R-squared: 0.9936
F-statistic: 1535 on 6 and 53 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,] 5.115014e-03 1.023003e-02 9.948850e-01
[2,] 9.489547e-04 1.897909e-03 9.990510e-01
[3,] 1.220856e-04 2.441713e-04 9.998779e-01
[4,] 2.053513e-05 4.107025e-05 9.999795e-01
[5,] 1.016850e-05 2.033701e-05 9.999898e-01
[6,] 9.896932e-06 1.979386e-05 9.999901e-01
[7,] 1.684063e-06 3.368126e-06 9.999983e-01
[8,] 2.779397e-07 5.558793e-07 9.999997e-01
[9,] 1.230490e-07 2.460979e-07 9.999999e-01
[10,] 1.420203e-07 2.840405e-07 9.999999e-01
[11,] 7.406184e-08 1.481237e-07 9.999999e-01
[12,] 7.365088e-08 1.473018e-07 9.999999e-01
[13,] 2.637966e-08 5.275931e-08 1.000000e+00
[14,] 6.788107e-09 1.357621e-08 1.000000e+00
[15,] 3.411102e-09 6.822205e-09 1.000000e+00
[16,] 7.404268e-10 1.480854e-09 1.000000e+00
[17,] 3.320227e-10 6.640454e-10 1.000000e+00
[18,] 1.084650e-10 2.169299e-10 1.000000e+00
[19,] 5.976534e-11 1.195307e-10 1.000000e+00
[20,] 1.177864e-11 2.355729e-11 1.000000e+00
[21,] 7.396146e-12 1.479229e-11 1.000000e+00
[22,] 9.765829e-12 1.953166e-11 1.000000e+00
[23,] 5.674769e-12 1.134954e-11 1.000000e+00
[24,] 1.369970e-12 2.739939e-12 1.000000e+00
[25,] 3.889829e-12 7.779658e-12 1.000000e+00
[26,] 7.031358e-10 1.406272e-09 1.000000e+00
[27,] 5.993342e-04 1.198668e-03 9.994007e-01
[28,] 5.833321e-04 1.166664e-03 9.994167e-01
[29,] 3.661352e-04 7.322705e-04 9.996339e-01
[30,] 3.363154e-04 6.726308e-04 9.996637e-01
[31,] 8.906809e-04 1.781362e-03 9.991093e-01
[32,] 2.115738e-03 4.231477e-03 9.978843e-01
[33,] 2.877732e-03 5.755464e-03 9.971223e-01
[34,] 3.868456e-03 7.736912e-03 9.961315e-01
[35,] 2.809814e-02 5.619628e-02 9.719019e-01
[36,] 5.142866e-02 1.028573e-01 9.485713e-01
[37,] 3.425182e-02 6.850364e-02 9.657482e-01
[38,] 2.205689e-02 4.411379e-02 9.779431e-01
[39,] 9.999675e-01 6.496241e-05 3.248121e-05
[40,] 9.999989e-01 2.290455e-06 1.145228e-06
[41,] 9.999713e-01 5.745716e-05 2.872858e-05
> postscript(file="/var/wessaorg/rcomp/tmp/1ou061321829715.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/2qlha1321829715.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/3douj1321829715.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/4kneo1321829715.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/5fudu1321829715.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
42.53991 205.35910 -101.87236 5810.11349 6748.93106 5805.09936
7 8 9 10 11 12
-506.20920 -3342.95542 -5392.77637 -2527.19246 -2340.80478 -2383.10577
13 14 15 16 17 18
-1479.94863 193.65567 689.28322 -1156.02050 -4425.91170 -3435.85707
19 20 21 22 23 24
408.19631 1716.06201 -711.05188 1646.36885 753.90570 -642.41077
25 26 27 28 29 30
2544.84078 3080.16681 1447.29238 1221.79451 -2064.21453 95.39259
31 32 33 34 35 36
267.43010 1083.64751 -418.69993 2473.31261 2798.67244 -6244.86928
37 38 39 40 41 42
-5010.75322 -1676.71584 -1558.85558 -969.80692 -1498.69995 881.41081
43 44 45 46 47 48
35.22304 2576.04884 3290.74103 3934.29425 4743.21324 3936.43846
49 50 51 52 53 54
-50.20427 245.09431 -815.21381 -1644.49886 -3350.25250 -1071.47420
55 56 57 58 59 60
-568.48988 3783.55846 -1782.34597 -1974.45488 -1085.98472 -2226.43559
> postscript(file="/var/wessaorg/rcomp/tmp/6psfw1321829715.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 42.53991 NA
1 205.35910 42.53991
2 -101.87236 205.35910
3 5810.11349 -101.87236
4 6748.93106 5810.11349
5 5805.09936 6748.93106
6 -506.20920 5805.09936
7 -3342.95542 -506.20920
8 -5392.77637 -3342.95542
9 -2527.19246 -5392.77637
10 -2340.80478 -2527.19246
11 -2383.10577 -2340.80478
12 -1479.94863 -2383.10577
13 193.65567 -1479.94863
14 689.28322 193.65567
15 -1156.02050 689.28322
16 -4425.91170 -1156.02050
17 -3435.85707 -4425.91170
18 408.19631 -3435.85707
19 1716.06201 408.19631
20 -711.05188 1716.06201
21 1646.36885 -711.05188
22 753.90570 1646.36885
23 -642.41077 753.90570
24 2544.84078 -642.41077
25 3080.16681 2544.84078
26 1447.29238 3080.16681
27 1221.79451 1447.29238
28 -2064.21453 1221.79451
29 95.39259 -2064.21453
30 267.43010 95.39259
31 1083.64751 267.43010
32 -418.69993 1083.64751
33 2473.31261 -418.69993
34 2798.67244 2473.31261
35 -6244.86928 2798.67244
36 -5010.75322 -6244.86928
37 -1676.71584 -5010.75322
38 -1558.85558 -1676.71584
39 -969.80692 -1558.85558
40 -1498.69995 -969.80692
41 881.41081 -1498.69995
42 35.22304 881.41081
43 2576.04884 35.22304
44 3290.74103 2576.04884
45 3934.29425 3290.74103
46 4743.21324 3934.29425
47 3936.43846 4743.21324
48 -50.20427 3936.43846
49 245.09431 -50.20427
50 -815.21381 245.09431
51 -1644.49886 -815.21381
52 -3350.25250 -1644.49886
53 -1071.47420 -3350.25250
54 -568.48988 -1071.47420
55 3783.55846 -568.48988
56 -1782.34597 3783.55846
57 -1974.45488 -1782.34597
58 -1085.98472 -1974.45488
59 -2226.43559 -1085.98472
60 NA -2226.43559
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 205.35910 42.53991
[2,] -101.87236 205.35910
[3,] 5810.11349 -101.87236
[4,] 6748.93106 5810.11349
[5,] 5805.09936 6748.93106
[6,] -506.20920 5805.09936
[7,] -3342.95542 -506.20920
[8,] -5392.77637 -3342.95542
[9,] -2527.19246 -5392.77637
[10,] -2340.80478 -2527.19246
[11,] -2383.10577 -2340.80478
[12,] -1479.94863 -2383.10577
[13,] 193.65567 -1479.94863
[14,] 689.28322 193.65567
[15,] -1156.02050 689.28322
[16,] -4425.91170 -1156.02050
[17,] -3435.85707 -4425.91170
[18,] 408.19631 -3435.85707
[19,] 1716.06201 408.19631
[20,] -711.05188 1716.06201
[21,] 1646.36885 -711.05188
[22,] 753.90570 1646.36885
[23,] -642.41077 753.90570
[24,] 2544.84078 -642.41077
[25,] 3080.16681 2544.84078
[26,] 1447.29238 3080.16681
[27,] 1221.79451 1447.29238
[28,] -2064.21453 1221.79451
[29,] 95.39259 -2064.21453
[30,] 267.43010 95.39259
[31,] 1083.64751 267.43010
[32,] -418.69993 1083.64751
[33,] 2473.31261 -418.69993
[34,] 2798.67244 2473.31261
[35,] -6244.86928 2798.67244
[36,] -5010.75322 -6244.86928
[37,] -1676.71584 -5010.75322
[38,] -1558.85558 -1676.71584
[39,] -969.80692 -1558.85558
[40,] -1498.69995 -969.80692
[41,] 881.41081 -1498.69995
[42,] 35.22304 881.41081
[43,] 2576.04884 35.22304
[44,] 3290.74103 2576.04884
[45,] 3934.29425 3290.74103
[46,] 4743.21324 3934.29425
[47,] 3936.43846 4743.21324
[48,] -50.20427 3936.43846
[49,] 245.09431 -50.20427
[50,] -815.21381 245.09431
[51,] -1644.49886 -815.21381
[52,] -3350.25250 -1644.49886
[53,] -1071.47420 -3350.25250
[54,] -568.48988 -1071.47420
[55,] 3783.55846 -568.48988
[56,] -1782.34597 3783.55846
[57,] -1974.45488 -1782.34597
[58,] -1085.98472 -1974.45488
[59,] -2226.43559 -1085.98472
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 205.35910 42.53991
2 -101.87236 205.35910
3 5810.11349 -101.87236
4 6748.93106 5810.11349
5 5805.09936 6748.93106
6 -506.20920 5805.09936
7 -3342.95542 -506.20920
8 -5392.77637 -3342.95542
9 -2527.19246 -5392.77637
10 -2340.80478 -2527.19246
11 -2383.10577 -2340.80478
12 -1479.94863 -2383.10577
13 193.65567 -1479.94863
14 689.28322 193.65567
15 -1156.02050 689.28322
16 -4425.91170 -1156.02050
17 -3435.85707 -4425.91170
18 408.19631 -3435.85707
19 1716.06201 408.19631
20 -711.05188 1716.06201
21 1646.36885 -711.05188
22 753.90570 1646.36885
23 -642.41077 753.90570
24 2544.84078 -642.41077
25 3080.16681 2544.84078
26 1447.29238 3080.16681
27 1221.79451 1447.29238
28 -2064.21453 1221.79451
29 95.39259 -2064.21453
30 267.43010 95.39259
31 1083.64751 267.43010
32 -418.69993 1083.64751
33 2473.31261 -418.69993
34 2798.67244 2473.31261
35 -6244.86928 2798.67244
36 -5010.75322 -6244.86928
37 -1676.71584 -5010.75322
38 -1558.85558 -1676.71584
39 -969.80692 -1558.85558
40 -1498.69995 -969.80692
41 881.41081 -1498.69995
42 35.22304 881.41081
43 2576.04884 35.22304
44 3290.74103 2576.04884
45 3934.29425 3290.74103
46 4743.21324 3934.29425
47 3936.43846 4743.21324
48 -50.20427 3936.43846
49 245.09431 -50.20427
50 -815.21381 245.09431
51 -1644.49886 -815.21381
52 -3350.25250 -1644.49886
53 -1071.47420 -3350.25250
54 -568.48988 -1071.47420
55 3783.55846 -568.48988
56 -1782.34597 3783.55846
57 -1974.45488 -1782.34597
58 -1085.98472 -1974.45488
59 -2226.43559 -1085.98472
> 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/7dldd1321829715.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/82ldr1321829715.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/9gex61321829715.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/10zfhd1321829715.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/113bxn1321829715.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/12zbce1321829715.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/136p9w1321829715.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/14iotw1321829715.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/15a6ad1321829715.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/16rn471321829715.tab")
+ }
>
> try(system("convert tmp/1ou061321829715.ps tmp/1ou061321829715.png",intern=TRUE))
character(0)
> try(system("convert tmp/2qlha1321829715.ps tmp/2qlha1321829715.png",intern=TRUE))
character(0)
> try(system("convert tmp/3douj1321829715.ps tmp/3douj1321829715.png",intern=TRUE))
character(0)
> try(system("convert tmp/4kneo1321829715.ps tmp/4kneo1321829715.png",intern=TRUE))
character(0)
> try(system("convert tmp/5fudu1321829715.ps tmp/5fudu1321829715.png",intern=TRUE))
character(0)
> try(system("convert tmp/6psfw1321829715.ps tmp/6psfw1321829715.png",intern=TRUE))
character(0)
> try(system("convert tmp/7dldd1321829715.ps tmp/7dldd1321829715.png",intern=TRUE))
character(0)
> try(system("convert tmp/82ldr1321829715.ps tmp/82ldr1321829715.png",intern=TRUE))
character(0)
> try(system("convert tmp/9gex61321829715.ps tmp/9gex61321829715.png",intern=TRUE))
character(0)
> try(system("convert tmp/10zfhd1321829715.ps tmp/10zfhd1321829715.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.346 0.515 3.906