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(631.923
+ ,9.911
+ ,58608
+ ,654.294
+ ,8.915
+ ,46865
+ ,671.833
+ ,9.452
+ ,51378
+ ,586.840
+ ,9.112
+ ,46235
+ ,600.969
+ ,8.472
+ ,47206
+ ,625.568
+ ,8.230
+ ,45382
+ ,558.110
+ ,8.384
+ ,41227
+ ,630.577
+ ,8.625
+ ,33795
+ ,628.654
+ ,8.221
+ ,31295
+ ,603.184
+ ,8.649
+ ,42625
+ ,656.255
+ ,8.625
+ ,33625
+ ,600.730
+ ,10.443
+ ,21538
+ ,670.326
+ ,10.357
+ ,56421
+ ,678.423
+ ,8.586
+ ,53152
+ ,641.502
+ ,8.892
+ ,53536
+ ,625.311
+ ,8.329
+ ,52408
+ ,628.177
+ ,8.101
+ ,41454
+ ,589.767
+ ,7.922
+ ,38271
+ ,582.471
+ ,8.120
+ ,35306
+ ,636.248
+ ,7.838
+ ,26414
+ ,599.885
+ ,7.735
+ ,31917
+ ,621.694
+ ,8.406
+ ,38030
+ ,637.406
+ ,8.209
+ ,27534
+ ,595.994
+ ,9.451
+ ,18387
+ ,696.308
+ ,10.041
+ ,50556
+ ,674.201
+ ,9.411
+ ,43901
+ ,648.861
+ ,10.405
+ ,48572
+ ,649.605
+ ,8.467
+ ,43899
+ ,672.392
+ ,8.464
+ ,37532
+ ,598.396
+ ,8.102
+ ,40357
+ ,613.177
+ ,7.627
+ ,35489
+ ,638.104
+ ,7.513
+ ,29027
+ ,615.632
+ ,7.510
+ ,34485
+ ,634.465
+ ,8.291
+ ,42598
+ ,638.686
+ ,8.064
+ ,30306
+ ,604.243
+ ,9.383
+ ,26451
+ ,706.669
+ ,9.706
+ ,47460
+ ,677.185
+ ,8.579
+ ,50104
+ ,644.328
+ ,9.474
+ ,61465
+ ,664.825
+ ,8.318
+ ,53726
+ ,605.707
+ ,8.213
+ ,39477
+ ,600.136
+ ,8.059
+ ,43895
+ ,612.166
+ ,9.111
+ ,31481
+ ,599.659
+ ,7.708
+ ,29896
+ ,634.210
+ ,7.680
+ ,33842
+ ,618.234
+ ,8.014
+ ,39120
+ ,613.576
+ ,8.007
+ ,33702
+ ,627.200
+ ,8.718
+ ,25094
+ ,668.973
+ ,9.486
+ ,51442
+ ,651.479
+ ,9.113
+ ,45594
+ ,619.661
+ ,9.025
+ ,52518
+ ,644.260
+ ,8.476
+ ,48564
+ ,579.936
+ ,7.952
+ ,41745
+ ,601.752
+ ,7.759
+ ,49585
+ ,595.376
+ ,7.835
+ ,32747
+ ,588.902
+ ,7.600
+ ,33379
+ ,634.341
+ ,7.651
+ ,35645
+ ,594.305
+ ,8.319
+ ,37034
+ ,606.200
+ ,8.812
+ ,35681
+ ,610.926
+ ,8.630
+ ,20972)
+ ,dim=c(3
+ ,60)
+ ,dimnames=list(c('WERKLOZEN'
+ ,'OVERLIJDENS'
+ ,'INSCHRIJVINGEN')
+ ,1:60))
> y <- array(NA,dim=c(3,60),dimnames=list(c('WERKLOZEN','OVERLIJDENS','INSCHRIJVINGEN'),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 = '1'
> 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
WERKLOZEN OVERLIJDENS INSCHRIJVINGEN
1 631.923 9.911 58608
2 654.294 8.915 46865
3 671.833 9.452 51378
4 586.840 9.112 46235
5 600.969 8.472 47206
6 625.568 8.230 45382
7 558.110 8.384 41227
8 630.577 8.625 33795
9 628.654 8.221 31295
10 603.184 8.649 42625
11 656.255 8.625 33625
12 600.730 10.443 21538
13 670.326 10.357 56421
14 678.423 8.586 53152
15 641.502 8.892 53536
16 625.311 8.329 52408
17 628.177 8.101 41454
18 589.767 7.922 38271
19 582.471 8.120 35306
20 636.248 7.838 26414
21 599.885 7.735 31917
22 621.694 8.406 38030
23 637.406 8.209 27534
24 595.994 9.451 18387
25 696.308 10.041 50556
26 674.201 9.411 43901
27 648.861 10.405 48572
28 649.605 8.467 43899
29 672.392 8.464 37532
30 598.396 8.102 40357
31 613.177 7.627 35489
32 638.104 7.513 29027
33 615.632 7.510 34485
34 634.465 8.291 42598
35 638.686 8.064 30306
36 604.243 9.383 26451
37 706.669 9.706 47460
38 677.185 8.579 50104
39 644.328 9.474 61465
40 664.825 8.318 53726
41 605.707 8.213 39477
42 600.136 8.059 43895
43 612.166 9.111 31481
44 599.659 7.708 29896
45 634.210 7.680 33842
46 618.234 8.014 39120
47 613.576 8.007 33702
48 627.200 8.718 25094
49 668.973 9.486 51442
50 651.479 9.113 45594
51 619.661 9.025 52518
52 644.260 8.476 48564
53 579.936 7.952 41745
54 601.752 7.759 49585
55 595.376 7.835 32747
56 588.902 7.600 33379
57 634.341 7.651 35645
58 594.305 8.319 37034
59 606.200 8.812 35681
60 610.926 8.630 20972
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) OVERLIJDENS INSCHRIJVINGEN
4.774e+02 1.265e+01 1.032e-03
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-67.903 -19.620 -0.415 20.298 57.500
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.774e+02 4.004e+01 11.924 < 2e-16 ***
OVERLIJDENS 1.265e+01 4.940e+00 2.561 0.01310 *
INSCHRIJVINGEN 1.032e-03 3.763e-04 2.741 0.00817 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 26.97 on 57 degrees of freedom
Multiple R-squared: 0.2713, Adjusted R-squared: 0.2458
F-statistic: 10.61 on 2 and 57 DF, p-value: 0.0001208
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.8937610 0.21247806 0.10623903
[2,] 0.9504550 0.09909003 0.04954501
[3,] 0.9499125 0.10017506 0.05008753
[4,] 0.9264030 0.14719407 0.07359704
[5,] 0.8986131 0.20277372 0.10138686
[6,] 0.8958814 0.20823723 0.10411861
[7,] 0.9435548 0.11289032 0.05644516
[8,] 0.9315841 0.13683183 0.06841592
[9,] 0.9581059 0.08378826 0.04189413
[10,] 0.9342490 0.13150203 0.06575101
[11,] 0.9058584 0.18828320 0.09414160
[12,] 0.8668892 0.26622167 0.13311084
[13,] 0.8567000 0.28659991 0.14329995
[14,] 0.8633768 0.27324645 0.13662322
[15,] 0.8931046 0.21379085 0.10689542
[16,] 0.8534456 0.29310878 0.14655439
[17,] 0.8035520 0.39289602 0.19644801
[18,] 0.8082706 0.38345878 0.19172939
[19,] 0.7809938 0.43801239 0.21900620
[20,] 0.8457244 0.30855112 0.15427556
[21,] 0.8603754 0.27924930 0.13962465
[22,] 0.8252754 0.34944922 0.17472461
[23,] 0.7990174 0.40196515 0.20098257
[24,] 0.8933658 0.21326841 0.10663420
[25,] 0.8826655 0.23466906 0.11733453
[26,] 0.8400476 0.31990479 0.15995240
[27,] 0.8772494 0.24550112 0.12275056
[28,] 0.8407055 0.31858904 0.15929452
[29,] 0.7933650 0.41327010 0.20663505
[30,] 0.8081625 0.38367499 0.19183749
[31,] 0.7960577 0.40788457 0.20394229
[32,] 0.9255972 0.14880555 0.07440278
[33,] 0.9662304 0.06753921 0.03376961
[34,] 0.9533438 0.09331250 0.04665625
[35,] 0.9712266 0.05754687 0.02877343
[36,] 0.9577017 0.08459653 0.04229826
[37,] 0.9470012 0.10599754 0.05299877
[38,] 0.9311883 0.13762346 0.06881173
[39,] 0.8940023 0.21199542 0.10599771
[40,] 0.9262415 0.14751707 0.07375854
[41,] 0.8909708 0.21805849 0.10902924
[42,] 0.8417175 0.31656505 0.15828252
[43,] 0.7841924 0.43161510 0.21580755
[44,] 0.7608548 0.47829036 0.23914518
[45,] 0.7559124 0.48817529 0.24408764
[46,] 0.6598447 0.68031052 0.34015526
[47,] 0.7901112 0.41977768 0.20988884
[48,] 0.7792283 0.44154337 0.22077169
[49,] 0.6261848 0.74763030 0.37381515
> postscript(file="/var/wessaorg/rcomp/tmp/1bkaf1324638463.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/22o5s1324638463.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/3kman1324638463.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/4ileo1324638463.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/54eu01324638463.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
-31.33904676 15.74662276 21.83621681 -53.54992227 -32.32535756 -2.78307146
7 8 9 10 11 12
-67.90339631 9.18093392 14.94814530 -27.62423786 35.03429627 -31.02352977
13 14 15 16 17 18
3.67720684 37.55277364 -3.63580949 -11.54022902 5.50992099 -27.35199306
19 20 21 22 23 24
-34.09453643 32.42277205 -8.31366851 -1.29989351 27.73160190 -19.95848745
25 26 27 28 29 30
39.70719724 32.43579225 -10.29849956 19.78521534 49.17800690 -23.15212791
31 32 33 34 35 36
2.66006030 35.69520543 7.63099788 8.21398126 27.98668430 -19.16752556
37 38 39 40 41 42
57.50022977 39.54748056 -16.35230113 26.75336746 -16.33772678 -24.51769781
43 44 45 46 47 48
-12.99188107 -6.11331923 24.72146325 -0.92474482 0.09471964 13.60276848
49 50 51 52 53 54
18.48003465 11.73764509 -26.10939641 9.51419922 -41.14613030 -24.97561967
55 56 57 58 59 60
-14.94403586 -19.09678307 23.35949481 -26.56076421 -19.50745533 2.69415389
> postscript(file="/var/wessaorg/rcomp/tmp/6cubz1324638463.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 -31.33904676 NA
1 15.74662276 -31.33904676
2 21.83621681 15.74662276
3 -53.54992227 21.83621681
4 -32.32535756 -53.54992227
5 -2.78307146 -32.32535756
6 -67.90339631 -2.78307146
7 9.18093392 -67.90339631
8 14.94814530 9.18093392
9 -27.62423786 14.94814530
10 35.03429627 -27.62423786
11 -31.02352977 35.03429627
12 3.67720684 -31.02352977
13 37.55277364 3.67720684
14 -3.63580949 37.55277364
15 -11.54022902 -3.63580949
16 5.50992099 -11.54022902
17 -27.35199306 5.50992099
18 -34.09453643 -27.35199306
19 32.42277205 -34.09453643
20 -8.31366851 32.42277205
21 -1.29989351 -8.31366851
22 27.73160190 -1.29989351
23 -19.95848745 27.73160190
24 39.70719724 -19.95848745
25 32.43579225 39.70719724
26 -10.29849956 32.43579225
27 19.78521534 -10.29849956
28 49.17800690 19.78521534
29 -23.15212791 49.17800690
30 2.66006030 -23.15212791
31 35.69520543 2.66006030
32 7.63099788 35.69520543
33 8.21398126 7.63099788
34 27.98668430 8.21398126
35 -19.16752556 27.98668430
36 57.50022977 -19.16752556
37 39.54748056 57.50022977
38 -16.35230113 39.54748056
39 26.75336746 -16.35230113
40 -16.33772678 26.75336746
41 -24.51769781 -16.33772678
42 -12.99188107 -24.51769781
43 -6.11331923 -12.99188107
44 24.72146325 -6.11331923
45 -0.92474482 24.72146325
46 0.09471964 -0.92474482
47 13.60276848 0.09471964
48 18.48003465 13.60276848
49 11.73764509 18.48003465
50 -26.10939641 11.73764509
51 9.51419922 -26.10939641
52 -41.14613030 9.51419922
53 -24.97561967 -41.14613030
54 -14.94403586 -24.97561967
55 -19.09678307 -14.94403586
56 23.35949481 -19.09678307
57 -26.56076421 23.35949481
58 -19.50745533 -26.56076421
59 2.69415389 -19.50745533
60 NA 2.69415389
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 15.74662276 -31.33904676
[2,] 21.83621681 15.74662276
[3,] -53.54992227 21.83621681
[4,] -32.32535756 -53.54992227
[5,] -2.78307146 -32.32535756
[6,] -67.90339631 -2.78307146
[7,] 9.18093392 -67.90339631
[8,] 14.94814530 9.18093392
[9,] -27.62423786 14.94814530
[10,] 35.03429627 -27.62423786
[11,] -31.02352977 35.03429627
[12,] 3.67720684 -31.02352977
[13,] 37.55277364 3.67720684
[14,] -3.63580949 37.55277364
[15,] -11.54022902 -3.63580949
[16,] 5.50992099 -11.54022902
[17,] -27.35199306 5.50992099
[18,] -34.09453643 -27.35199306
[19,] 32.42277205 -34.09453643
[20,] -8.31366851 32.42277205
[21,] -1.29989351 -8.31366851
[22,] 27.73160190 -1.29989351
[23,] -19.95848745 27.73160190
[24,] 39.70719724 -19.95848745
[25,] 32.43579225 39.70719724
[26,] -10.29849956 32.43579225
[27,] 19.78521534 -10.29849956
[28,] 49.17800690 19.78521534
[29,] -23.15212791 49.17800690
[30,] 2.66006030 -23.15212791
[31,] 35.69520543 2.66006030
[32,] 7.63099788 35.69520543
[33,] 8.21398126 7.63099788
[34,] 27.98668430 8.21398126
[35,] -19.16752556 27.98668430
[36,] 57.50022977 -19.16752556
[37,] 39.54748056 57.50022977
[38,] -16.35230113 39.54748056
[39,] 26.75336746 -16.35230113
[40,] -16.33772678 26.75336746
[41,] -24.51769781 -16.33772678
[42,] -12.99188107 -24.51769781
[43,] -6.11331923 -12.99188107
[44,] 24.72146325 -6.11331923
[45,] -0.92474482 24.72146325
[46,] 0.09471964 -0.92474482
[47,] 13.60276848 0.09471964
[48,] 18.48003465 13.60276848
[49,] 11.73764509 18.48003465
[50,] -26.10939641 11.73764509
[51,] 9.51419922 -26.10939641
[52,] -41.14613030 9.51419922
[53,] -24.97561967 -41.14613030
[54,] -14.94403586 -24.97561967
[55,] -19.09678307 -14.94403586
[56,] 23.35949481 -19.09678307
[57,] -26.56076421 23.35949481
[58,] -19.50745533 -26.56076421
[59,] 2.69415389 -19.50745533
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 15.74662276 -31.33904676
2 21.83621681 15.74662276
3 -53.54992227 21.83621681
4 -32.32535756 -53.54992227
5 -2.78307146 -32.32535756
6 -67.90339631 -2.78307146
7 9.18093392 -67.90339631
8 14.94814530 9.18093392
9 -27.62423786 14.94814530
10 35.03429627 -27.62423786
11 -31.02352977 35.03429627
12 3.67720684 -31.02352977
13 37.55277364 3.67720684
14 -3.63580949 37.55277364
15 -11.54022902 -3.63580949
16 5.50992099 -11.54022902
17 -27.35199306 5.50992099
18 -34.09453643 -27.35199306
19 32.42277205 -34.09453643
20 -8.31366851 32.42277205
21 -1.29989351 -8.31366851
22 27.73160190 -1.29989351
23 -19.95848745 27.73160190
24 39.70719724 -19.95848745
25 32.43579225 39.70719724
26 -10.29849956 32.43579225
27 19.78521534 -10.29849956
28 49.17800690 19.78521534
29 -23.15212791 49.17800690
30 2.66006030 -23.15212791
31 35.69520543 2.66006030
32 7.63099788 35.69520543
33 8.21398126 7.63099788
34 27.98668430 8.21398126
35 -19.16752556 27.98668430
36 57.50022977 -19.16752556
37 39.54748056 57.50022977
38 -16.35230113 39.54748056
39 26.75336746 -16.35230113
40 -16.33772678 26.75336746
41 -24.51769781 -16.33772678
42 -12.99188107 -24.51769781
43 -6.11331923 -12.99188107
44 24.72146325 -6.11331923
45 -0.92474482 24.72146325
46 0.09471964 -0.92474482
47 13.60276848 0.09471964
48 18.48003465 13.60276848
49 11.73764509 18.48003465
50 -26.10939641 11.73764509
51 9.51419922 -26.10939641
52 -41.14613030 9.51419922
53 -24.97561967 -41.14613030
54 -14.94403586 -24.97561967
55 -19.09678307 -14.94403586
56 23.35949481 -19.09678307
57 -26.56076421 23.35949481
58 -19.50745533 -26.56076421
59 2.69415389 -19.50745533
> 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/77tzl1324638463.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/8g0pr1324638463.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/9lz7o1324638463.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/10ocg81324638463.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/11zrwi1324638463.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/12da6h1324638463.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/13o1h01324638463.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/14vs031324638463.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/15m0p31324638463.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/161eum1324638463.tab")
+ }
>
> try(system("convert tmp/1bkaf1324638463.ps tmp/1bkaf1324638463.png",intern=TRUE))
character(0)
> try(system("convert tmp/22o5s1324638463.ps tmp/22o5s1324638463.png",intern=TRUE))
character(0)
> try(system("convert tmp/3kman1324638463.ps tmp/3kman1324638463.png",intern=TRUE))
character(0)
> try(system("convert tmp/4ileo1324638463.ps tmp/4ileo1324638463.png",intern=TRUE))
character(0)
> try(system("convert tmp/54eu01324638463.ps tmp/54eu01324638463.png",intern=TRUE))
character(0)
> try(system("convert tmp/6cubz1324638463.ps tmp/6cubz1324638463.png",intern=TRUE))
character(0)
> try(system("convert tmp/77tzl1324638463.ps tmp/77tzl1324638463.png",intern=TRUE))
character(0)
> try(system("convert tmp/8g0pr1324638463.ps tmp/8g0pr1324638463.png",intern=TRUE))
character(0)
> try(system("convert tmp/9lz7o1324638463.ps tmp/9lz7o1324638463.png",intern=TRUE))
character(0)
> try(system("convert tmp/10ocg81324638463.ps tmp/10ocg81324638463.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.108 0.535 3.673