R version 2.12.0 (2010-10-15)
Copyright (C) 2010 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(349774
+ ,100.7
+ ,98.7
+ ,351112
+ ,99.7937
+ ,99.2922
+ ,352582
+ ,100.2926685
+ ,99.8879532
+ ,354211
+ ,102.0979365
+ ,102.9844797
+ ,355695
+ ,102.8126221
+ ,106.5889365
+ ,358062
+ ,103.2238726
+ ,111.3854387
+ ,360249
+ ,104.0496636
+ ,112.8334494
+ ,362373
+ ,102.4889186
+ ,113.0591163
+ ,364607
+ ,102.7963854
+ ,111.4762887
+ ,367170
+ ,105.1607022
+ ,116.8271505
+ ,369253
+ ,107.4742377
+ ,120.7992736
+ ,372173
+ ,113.062898
+ ,122.732062
+ ,374406
+ ,121.3164896
+ ,125.5548994
+ ,376728
+ ,126.5330986
+ ,128.6937719
+ ,378832
+ ,132.3536212
+ ,135.7719294
+ ,381109
+ ,132.0889139
+ ,139.9808592
+ ,383366
+ ,139.2217153
+ ,145.4401127
+ ,385808
+ ,136.437281
+ ,145.876433
+ ,388012
+ ,136.7101555
+ ,150.6903553
+ ,390714
+ ,144.3659243
+ ,161.3893706
+ ,393057
+ ,160.8236396
+ ,158.96853
+ ,395384
+ ,157.9288141
+ ,173.9115718
+ ,396134
+ ,165.3514684
+ ,174.0854834
+ ,396948
+ ,164.6900625
+ ,181.0489027
+ ,397499
+ ,173.5833259
+ ,178.695267
+ ,396441
+ ,170.979576
+ ,191.9187167
+ ,395202
+ ,171.3215351
+ ,196.7166847
+ ,392817
+ ,179.0310042
+ ,208.1262524
+ ,391455
+ ,184.0438723
+ ,226.4413626
+ ,391036
+ ,185.1481356
+ ,227.5735694
+ ,390741
+ ,185.1481356
+ ,236.221365
+ ,390853
+ ,184.0372468
+ ,235.7489223
+ ,390270
+ ,188.4541407
+ ,244.4716324
+ ,389527
+ ,202.9651095
+ ,257.9175722
+ ,389611
+ ,208.0392372
+ ,282.4197416
+ ,390544
+ ,211.7839435
+ ,286.6560377
+ ,391684
+ ,200.7711785
+ ,288.3759739
+ ,393854
+ ,205.9912291
+ ,297.8923811
+ ,394869
+ ,213.6129046
+ ,299.9776277
+ ,396623
+ ,211.2631626
+ ,301.4775159
+ ,397981
+ ,214.0095837
+ ,314.1395715
+ ,400169
+ ,213.7955741
+ ,311.626455
+ ,402390
+ ,213.1541874
+ ,321.2868751
+ ,403789
+ ,210.383183
+ ,320.6443013
+ ,406203
+ ,216.9050617
+ ,328.6604088
+ ,407742
+ ,213.0007706
+ ,329.9750505
+ ,409045
+ ,207.4627505
+ ,322.7155994
+ ,410108
+ ,231.5284296
+ ,331.4289206
+ ,411676
+ ,231.5284296
+ ,331.0974916
+ ,412786
+ ,223.8879914
+ ,342.0237089
+ ,412931
+ ,224.3357674
+ ,358.0988232
+ ,413654
+ ,227.9251397
+ ,365.6188985
+ ,413750
+ ,231.1160916
+ ,356.8440449
+ ,412324
+ ,216.5557778
+ ,364.6946139
+ ,410074
+ ,210.7087718
+ ,362.1417516
+ ,405189
+ ,237.0473683
+ ,349.828932
+ ,398441
+ ,236.0991789
+ ,354.3767081
+ ,392869
+ ,247.1958403
+ ,382.7268448
+ ,389881
+ ,243.9822943
+ ,407.6040897
+ ,388275
+ ,240.3225599
+ ,430.0223146
+ ,387965
+ ,240.803205
+ ,449.8033411
+ ,389869
+ ,248.2681044
+ ,455.2009812
+ ,389649
+ ,249.2611768
+ ,464.7602018
+ ,390383
+ ,244.0266921
+ ,474.9849263
+ ,391648
+ ,244.7587722
+ ,472.1350167
+ ,393048
+ ,241.5769081
+ ,471.6628817
+ ,393888
+ ,235.7790623
+ ,486.284431)
+ ,dim=c(3
+ ,67)
+ ,dimnames=list(c('Werkgelegenheid'
+ ,'Uurloon'
+ ,'Productiviteit')
+ ,1:67))
> y <- array(NA,dim=c(3,67),dimnames=list(c('Werkgelegenheid','Uurloon','Productiviteit'),1:67))
> 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 = '3'
> 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
Productiviteit Werkgelegenheid Uurloon
1 98.70000 349774 100.7000
2 99.29220 351112 99.7937
3 99.88795 352582 100.2927
4 102.98448 354211 102.0979
5 106.58894 355695 102.8126
6 111.38544 358062 103.2239
7 112.83345 360249 104.0497
8 113.05912 362373 102.4889
9 111.47629 364607 102.7964
10 116.82715 367170 105.1607
11 120.79927 369253 107.4742
12 122.73206 372173 113.0629
13 125.55490 374406 121.3165
14 128.69377 376728 126.5331
15 135.77193 378832 132.3536
16 139.98086 381109 132.0889
17 145.44011 383366 139.2217
18 145.87643 385808 136.4373
19 150.69036 388012 136.7102
20 161.38937 390714 144.3659
21 158.96853 393057 160.8236
22 173.91157 395384 157.9288
23 174.08548 396134 165.3515
24 181.04890 396948 164.6901
25 178.69527 397499 173.5833
26 191.91872 396441 170.9796
27 196.71668 395202 171.3215
28 208.12625 392817 179.0310
29 226.44136 391455 184.0439
30 227.57357 391036 185.1481
31 236.22136 390741 185.1481
32 235.74892 390853 184.0372
33 244.47163 390270 188.4541
34 257.91757 389527 202.9651
35 282.41974 389611 208.0392
36 286.65604 390544 211.7839
37 288.37597 391684 200.7712
38 297.89238 393854 205.9912
39 299.97763 394869 213.6129
40 301.47752 396623 211.2632
41 314.13957 397981 214.0096
42 311.62646 400169 213.7956
43 321.28688 402390 213.1542
44 320.64430 403789 210.3832
45 328.66041 406203 216.9051
46 329.97505 407742 213.0008
47 322.71560 409045 207.4628
48 331.42892 410108 231.5284
49 331.09749 411676 231.5284
50 342.02371 412786 223.8880
51 358.09882 412931 224.3358
52 365.61890 413654 227.9251
53 356.84404 413750 231.1161
54 364.69461 412324 216.5558
55 362.14175 410074 210.7088
56 349.82893 405189 237.0474
57 354.37671 398441 236.0992
58 382.72684 392869 247.1958
59 407.60409 389881 243.9823
60 430.02231 388275 240.3226
61 449.80334 387965 240.8032
62 455.20098 389869 248.2681
63 464.76020 389649 249.2612
64 474.98493 390383 244.0267
65 472.13502 391648 244.7588
66 471.66288 393048 241.5769
67 486.28443 393888 235.7791
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Werkgelegenheid Uurloon
538.031863 -0.002027 2.793732
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-57.442 -22.341 -2.176 18.834 88.093
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.380e+02 1.286e+02 4.185 8.88e-05 ***
Werkgelegenheid -2.027e-03 3.753e-04 -5.402 1.03e-06 ***
Uurloon 2.794e+00 1.249e-01 22.372 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 30.38 on 64 degrees of freedom
Multiple R-squared: 0.9362, Adjusted R-squared: 0.9342
F-statistic: 469.4 on 2 and 64 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.213294e-04 2.426589e-04 0.99987867
[2,] 4.664735e-06 9.329470e-06 0.99999534
[3,] 1.651173e-07 3.302346e-07 0.99999983
[4,] 6.971808e-08 1.394362e-07 0.99999993
[5,] 8.347231e-09 1.669446e-08 0.99999999
[6,] 8.978571e-10 1.795714e-09 1.00000000
[7,] 4.363271e-10 8.726541e-10 1.00000000
[8,] 3.902982e-11 7.805965e-11 1.00000000
[9,] 2.486930e-12 4.973860e-12 1.00000000
[10,] 5.088025e-13 1.017605e-12 1.00000000
[11,] 1.439167e-13 2.878335e-13 1.00000000
[12,] 3.647585e-14 7.295169e-14 1.00000000
[13,] 4.324279e-15 8.648557e-15 1.00000000
[14,] 1.924508e-15 3.849016e-15 1.00000000
[15,] 2.095486e-14 4.190972e-14 1.00000000
[16,] 3.356813e-15 6.713627e-15 1.00000000
[17,] 3.210425e-14 6.420850e-14 1.00000000
[18,] 7.277150e-15 1.455430e-14 1.00000000
[19,] 1.344857e-14 2.689715e-14 1.00000000
[20,] 1.752985e-15 3.505971e-15 1.00000000
[21,] 3.129936e-14 6.259871e-14 1.00000000
[22,] 4.233665e-13 8.467329e-13 1.00000000
[23,] 2.873918e-12 5.747835e-12 1.00000000
[24,] 5.549068e-11 1.109814e-10 1.00000000
[25,] 7.159826e-11 1.431965e-10 1.00000000
[26,] 1.759807e-10 3.519614e-10 1.00000000
[27,] 2.204614e-10 4.409228e-10 1.00000000
[28,] 1.761502e-10 3.523005e-10 1.00000000
[29,] 7.178129e-11 1.435626e-10 1.00000000
[30,] 7.559910e-11 1.511982e-10 1.00000000
[31,] 7.325367e-11 1.465073e-10 1.00000000
[32,] 8.260268e-10 1.652054e-09 1.00000000
[33,] 3.982808e-09 7.965615e-09 1.00000000
[34,] 5.658033e-09 1.131607e-08 0.99999999
[35,] 1.411504e-08 2.823008e-08 0.99999999
[36,] 6.870027e-08 1.374005e-07 0.99999993
[37,] 3.113714e-07 6.227428e-07 0.99999969
[38,] 2.818062e-06 5.636124e-06 0.99999718
[39,] 2.622217e-05 5.244435e-05 0.99997378
[40,] 8.485417e-05 1.697083e-04 0.99991515
[41,] 3.557794e-04 7.115587e-04 0.99964422
[42,] 2.138415e-03 4.276829e-03 0.99786159
[43,] 1.541387e-03 3.082773e-03 0.99845861
[44,] 1.009718e-03 2.019437e-03 0.99899028
[45,] 9.722241e-04 1.944448e-03 0.99902778
[46,] 1.352679e-03 2.705358e-03 0.99864732
[47,] 1.918292e-03 3.836585e-03 0.99808171
[48,] 1.956343e-03 3.912686e-03 0.99804366
[49,] 4.559732e-03 9.119465e-03 0.99544027
[50,] 7.494989e-03 1.498998e-02 0.99250501
[51,] 3.819345e-03 7.638691e-03 0.99618065
[52,] 3.463867e-02 6.927734e-02 0.96536133
[53,] 5.009997e-01 9.980007e-01 0.49900035
[54,] 9.556777e-01 8.864469e-02 0.04432235
[55,] 9.814102e-01 3.717951e-02 0.01858976
[56,] 9.838217e-01 3.235656e-02 0.01617828
> postscript(file="/var/www/rcomp/tmp/1blvl1322144170.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/www/rcomp/tmp/2g70n1322144170.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/www/rcomp/tmp/3lsw31322144170.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/www/rcomp/tmp/4yxtk1322144170.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/www/rcomp/tmp/5bk0a1322144170.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 = 67
Frequency = 1
1 2 3 4 5 6
-11.5507591 -5.7140221 -3.5320672 -2.1764432 2.4399423 10.8862285
7 8 9 10 11 12
14.4609873 23.3530224 25.4402867 29.3819466 31.1136152 23.3530092
13 14 15 16 17 18
7.6445673 0.9171092 -4.0001962 5.5645023 -4.3276799 8.8383619
19 20 21 22 23 24
17.3581977 12.1469116 -31.5023250 -3.7543020 -22.7967950 -12.3353321
25 26 27 28 29 30
-38.4173004 -20.0645940 -18.7338392 -33.6976627 -32.1483967 -34.9506601
31 32 33 34 35 36
-26.9009291 -24.0427843 -28.8416302 -57.4417616 -46.9450494 -51.2789537
37 38 39 40 41 42
-16.4811376 -17.1488311 -34.2987566 -22.6783688 -14.9359538 -12.4153705
43 44 45 46 47 48
3.5396287 13.4747440 8.1644621 23.5067196 34.3606337 -22.0040458
49 50 51 52 53 54
-19.1566097 15.3652897 30.4834017 29.4414973 11.9466033 57.5838054
55 56 57 58 59 60
66.8044022 -28.9949436 -35.4786534 -49.4259402 -21.6285929 7.7580456
61 62 63 64 65 66
25.5678038 13.9705652 20.3093934 46.6459315 44.3153681 55.5707814
67
88.0929222
> postscript(file="/var/www/rcomp/tmp/69ss41322144170.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 = 67
Frequency = 1
lag(myerror, k = 1) myerror
0 -11.5507591 NA
1 -5.7140221 -11.5507591
2 -3.5320672 -5.7140221
3 -2.1764432 -3.5320672
4 2.4399423 -2.1764432
5 10.8862285 2.4399423
6 14.4609873 10.8862285
7 23.3530224 14.4609873
8 25.4402867 23.3530224
9 29.3819466 25.4402867
10 31.1136152 29.3819466
11 23.3530092 31.1136152
12 7.6445673 23.3530092
13 0.9171092 7.6445673
14 -4.0001962 0.9171092
15 5.5645023 -4.0001962
16 -4.3276799 5.5645023
17 8.8383619 -4.3276799
18 17.3581977 8.8383619
19 12.1469116 17.3581977
20 -31.5023250 12.1469116
21 -3.7543020 -31.5023250
22 -22.7967950 -3.7543020
23 -12.3353321 -22.7967950
24 -38.4173004 -12.3353321
25 -20.0645940 -38.4173004
26 -18.7338392 -20.0645940
27 -33.6976627 -18.7338392
28 -32.1483967 -33.6976627
29 -34.9506601 -32.1483967
30 -26.9009291 -34.9506601
31 -24.0427843 -26.9009291
32 -28.8416302 -24.0427843
33 -57.4417616 -28.8416302
34 -46.9450494 -57.4417616
35 -51.2789537 -46.9450494
36 -16.4811376 -51.2789537
37 -17.1488311 -16.4811376
38 -34.2987566 -17.1488311
39 -22.6783688 -34.2987566
40 -14.9359538 -22.6783688
41 -12.4153705 -14.9359538
42 3.5396287 -12.4153705
43 13.4747440 3.5396287
44 8.1644621 13.4747440
45 23.5067196 8.1644621
46 34.3606337 23.5067196
47 -22.0040458 34.3606337
48 -19.1566097 -22.0040458
49 15.3652897 -19.1566097
50 30.4834017 15.3652897
51 29.4414973 30.4834017
52 11.9466033 29.4414973
53 57.5838054 11.9466033
54 66.8044022 57.5838054
55 -28.9949436 66.8044022
56 -35.4786534 -28.9949436
57 -49.4259402 -35.4786534
58 -21.6285929 -49.4259402
59 7.7580456 -21.6285929
60 25.5678038 7.7580456
61 13.9705652 25.5678038
62 20.3093934 13.9705652
63 46.6459315 20.3093934
64 44.3153681 46.6459315
65 55.5707814 44.3153681
66 88.0929222 55.5707814
67 NA 88.0929222
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -5.7140221 -11.5507591
[2,] -3.5320672 -5.7140221
[3,] -2.1764432 -3.5320672
[4,] 2.4399423 -2.1764432
[5,] 10.8862285 2.4399423
[6,] 14.4609873 10.8862285
[7,] 23.3530224 14.4609873
[8,] 25.4402867 23.3530224
[9,] 29.3819466 25.4402867
[10,] 31.1136152 29.3819466
[11,] 23.3530092 31.1136152
[12,] 7.6445673 23.3530092
[13,] 0.9171092 7.6445673
[14,] -4.0001962 0.9171092
[15,] 5.5645023 -4.0001962
[16,] -4.3276799 5.5645023
[17,] 8.8383619 -4.3276799
[18,] 17.3581977 8.8383619
[19,] 12.1469116 17.3581977
[20,] -31.5023250 12.1469116
[21,] -3.7543020 -31.5023250
[22,] -22.7967950 -3.7543020
[23,] -12.3353321 -22.7967950
[24,] -38.4173004 -12.3353321
[25,] -20.0645940 -38.4173004
[26,] -18.7338392 -20.0645940
[27,] -33.6976627 -18.7338392
[28,] -32.1483967 -33.6976627
[29,] -34.9506601 -32.1483967
[30,] -26.9009291 -34.9506601
[31,] -24.0427843 -26.9009291
[32,] -28.8416302 -24.0427843
[33,] -57.4417616 -28.8416302
[34,] -46.9450494 -57.4417616
[35,] -51.2789537 -46.9450494
[36,] -16.4811376 -51.2789537
[37,] -17.1488311 -16.4811376
[38,] -34.2987566 -17.1488311
[39,] -22.6783688 -34.2987566
[40,] -14.9359538 -22.6783688
[41,] -12.4153705 -14.9359538
[42,] 3.5396287 -12.4153705
[43,] 13.4747440 3.5396287
[44,] 8.1644621 13.4747440
[45,] 23.5067196 8.1644621
[46,] 34.3606337 23.5067196
[47,] -22.0040458 34.3606337
[48,] -19.1566097 -22.0040458
[49,] 15.3652897 -19.1566097
[50,] 30.4834017 15.3652897
[51,] 29.4414973 30.4834017
[52,] 11.9466033 29.4414973
[53,] 57.5838054 11.9466033
[54,] 66.8044022 57.5838054
[55,] -28.9949436 66.8044022
[56,] -35.4786534 -28.9949436
[57,] -49.4259402 -35.4786534
[58,] -21.6285929 -49.4259402
[59,] 7.7580456 -21.6285929
[60,] 25.5678038 7.7580456
[61,] 13.9705652 25.5678038
[62,] 20.3093934 13.9705652
[63,] 46.6459315 20.3093934
[64,] 44.3153681 46.6459315
[65,] 55.5707814 44.3153681
[66,] 88.0929222 55.5707814
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -5.7140221 -11.5507591
2 -3.5320672 -5.7140221
3 -2.1764432 -3.5320672
4 2.4399423 -2.1764432
5 10.8862285 2.4399423
6 14.4609873 10.8862285
7 23.3530224 14.4609873
8 25.4402867 23.3530224
9 29.3819466 25.4402867
10 31.1136152 29.3819466
11 23.3530092 31.1136152
12 7.6445673 23.3530092
13 0.9171092 7.6445673
14 -4.0001962 0.9171092
15 5.5645023 -4.0001962
16 -4.3276799 5.5645023
17 8.8383619 -4.3276799
18 17.3581977 8.8383619
19 12.1469116 17.3581977
20 -31.5023250 12.1469116
21 -3.7543020 -31.5023250
22 -22.7967950 -3.7543020
23 -12.3353321 -22.7967950
24 -38.4173004 -12.3353321
25 -20.0645940 -38.4173004
26 -18.7338392 -20.0645940
27 -33.6976627 -18.7338392
28 -32.1483967 -33.6976627
29 -34.9506601 -32.1483967
30 -26.9009291 -34.9506601
31 -24.0427843 -26.9009291
32 -28.8416302 -24.0427843
33 -57.4417616 -28.8416302
34 -46.9450494 -57.4417616
35 -51.2789537 -46.9450494
36 -16.4811376 -51.2789537
37 -17.1488311 -16.4811376
38 -34.2987566 -17.1488311
39 -22.6783688 -34.2987566
40 -14.9359538 -22.6783688
41 -12.4153705 -14.9359538
42 3.5396287 -12.4153705
43 13.4747440 3.5396287
44 8.1644621 13.4747440
45 23.5067196 8.1644621
46 34.3606337 23.5067196
47 -22.0040458 34.3606337
48 -19.1566097 -22.0040458
49 15.3652897 -19.1566097
50 30.4834017 15.3652897
51 29.4414973 30.4834017
52 11.9466033 29.4414973
53 57.5838054 11.9466033
54 66.8044022 57.5838054
55 -28.9949436 66.8044022
56 -35.4786534 -28.9949436
57 -49.4259402 -35.4786534
58 -21.6285929 -49.4259402
59 7.7580456 -21.6285929
60 25.5678038 7.7580456
61 13.9705652 25.5678038
62 20.3093934 13.9705652
63 46.6459315 20.3093934
64 44.3153681 46.6459315
65 55.5707814 44.3153681
66 88.0929222 55.5707814
> 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/www/rcomp/tmp/7n3bs1322144170.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/www/rcomp/tmp/896sl1322144170.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/www/rcomp/tmp/9sn681322144170.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/www/rcomp/tmp/10u6wg1322144170.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/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/www/rcomp/tmp/11zsky1322144170.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/www/rcomp/tmp/12qdi11322144170.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/www/rcomp/tmp/13xjc61322144170.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/www/rcomp/tmp/14c4az1322144170.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/www/rcomp/tmp/158msh1322144170.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/www/rcomp/tmp/16sk2k1322144170.tab")
+ }
>
> try(system("convert tmp/1blvl1322144170.ps tmp/1blvl1322144170.png",intern=TRUE))
character(0)
> try(system("convert tmp/2g70n1322144170.ps tmp/2g70n1322144170.png",intern=TRUE))
character(0)
> try(system("convert tmp/3lsw31322144170.ps tmp/3lsw31322144170.png",intern=TRUE))
character(0)
> try(system("convert tmp/4yxtk1322144170.ps tmp/4yxtk1322144170.png",intern=TRUE))
character(0)
> try(system("convert tmp/5bk0a1322144170.ps tmp/5bk0a1322144170.png",intern=TRUE))
character(0)
> try(system("convert tmp/69ss41322144170.ps tmp/69ss41322144170.png",intern=TRUE))
character(0)
> try(system("convert tmp/7n3bs1322144170.ps tmp/7n3bs1322144170.png",intern=TRUE))
character(0)
> try(system("convert tmp/896sl1322144170.ps tmp/896sl1322144170.png",intern=TRUE))
character(0)
> try(system("convert tmp/9sn681322144170.ps tmp/9sn681322144170.png",intern=TRUE))
character(0)
> try(system("convert tmp/10u6wg1322144170.ps tmp/10u6wg1322144170.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.30 0.25 4.55