R version 2.12.1 (2010-12-16)
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(12008.00
+ ,4.00
+ ,9169.00
+ ,5.90
+ ,8788.00
+ ,7.10
+ ,8417.00
+ ,10.50
+ ,8247.00
+ ,15.10
+ ,8197.00
+ ,16.80
+ ,8236.00
+ ,15.30
+ ,8253.00
+ ,18.40
+ ,7733.00
+ ,16.10
+ ,8366.00
+ ,11.30
+ ,8626.00
+ ,7.90
+ ,8863.00
+ ,5.60
+ ,10102.00
+ ,3.40
+ ,8463.00
+ ,4.80
+ ,9114.00
+ ,6.50
+ ,8563.00
+ ,8.50
+ ,8872.00
+ ,15.10
+ ,8301.00
+ ,15.70
+ ,8301.00
+ ,18.70
+ ,8278.00
+ ,19.20
+ ,7736.00
+ ,12.90
+ ,7973.00
+ ,14.40
+ ,8268.00
+ ,6.20
+ ,9476.00
+ ,3.30
+ ,11100.00
+ ,4.60
+ ,8962.00
+ ,7.10
+ ,9173.00
+ ,7.80
+ ,8738.00
+ ,9.90
+ ,8459.00
+ ,13.60
+ ,8078.00
+ ,17.10
+ ,8411.00
+ ,17.80
+ ,8291.00
+ ,18.60
+ ,7810.00
+ ,14.70
+ ,8616.00
+ ,10.50
+ ,8312.00
+ ,8.60
+ ,9692.00
+ ,4.40
+ ,9911.00
+ ,2.30
+ ,8915.00
+ ,2.80
+ ,9452.00
+ ,8.80
+ ,9112.00
+ ,10.70
+ ,8472.00
+ ,13.90
+ ,8230.00
+ ,19.30
+ ,8384.00
+ ,19.50
+ ,8625.00
+ ,20.40
+ ,8221.00
+ ,15.30
+ ,8649.00
+ ,7.90
+ ,8625.00
+ ,8.30
+ ,10443.00
+ ,4.50
+ ,10357.00
+ ,3.20
+ ,8586.00
+ ,5.00
+ ,8892.00
+ ,6.60
+ ,8329.00
+ ,11.10
+ ,8101.00
+ ,12.80
+ ,7922.00
+ ,16.30
+ ,8120.00
+ ,17.40
+ ,7838.00
+ ,18.90
+ ,7735.00
+ ,15.80
+ ,8406.00
+ ,11.70
+ ,8209.00
+ ,6.40
+ ,9451.00
+ ,2.90
+ ,10041.00
+ ,4.70
+ ,9411.00
+ ,2.40
+ ,10405.00
+ ,7.20
+ ,8467.00
+ ,10.70
+ ,8464.00
+ ,13.40
+ ,8102.00
+ ,18.30
+ ,7627.00
+ ,18.40
+ ,7513.00
+ ,16.80
+ ,7510.00
+ ,16.60
+ ,8291.00
+ ,14.10
+ ,8064.00
+ ,6.10
+ ,9383.00
+ ,3.50
+ ,9706.00
+ ,1.70
+ ,8579.00
+ ,2.30
+ ,9474.00
+ ,4.50
+ ,8318.00
+ ,9.30
+ ,8213.00
+ ,14.20
+ ,8059.00
+ ,17.30
+ ,9111.00
+ ,23.00
+ ,7708.00
+ ,16.30
+ ,7680.00
+ ,18.40
+ ,8014.00
+ ,14.20
+ ,8007.00
+ ,9.10
+ ,8718.00
+ ,5.90
+ ,9486.00
+ ,7.20
+ ,9113.00
+ ,6.80
+ ,9025.00
+ ,8.00
+ ,8476.00
+ ,14.30
+ ,7952.00
+ ,14.60
+ ,7759.00
+ ,17.50
+ ,7835.00
+ ,17.20
+ ,7600.00
+ ,17.20
+ ,7651.00
+ ,14.10
+ ,8319.00
+ ,10.40
+ ,8812.00
+ ,6.80
+ ,8630.00
+ ,4.10)
+ ,dim=c(2
+ ,96)
+ ,dimnames=list(c('Sterftecijfers'
+ ,'Temperatuur')
+ ,1:96))
> y <- array(NA,dim=c(2,96),dimnames=list(c('Sterftecijfers','Temperatuur'),1:96))
> 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
Sterftecijfers Temperatuur
1 12008 4.0
2 9169 5.9
3 8788 7.1
4 8417 10.5
5 8247 15.1
6 8197 16.8
7 8236 15.3
8 8253 18.4
9 7733 16.1
10 8366 11.3
11 8626 7.9
12 8863 5.6
13 10102 3.4
14 8463 4.8
15 9114 6.5
16 8563 8.5
17 8872 15.1
18 8301 15.7
19 8301 18.7
20 8278 19.2
21 7736 12.9
22 7973 14.4
23 8268 6.2
24 9476 3.3
25 11100 4.6
26 8962 7.1
27 9173 7.8
28 8738 9.9
29 8459 13.6
30 8078 17.1
31 8411 17.8
32 8291 18.6
33 7810 14.7
34 8616 10.5
35 8312 8.6
36 9692 4.4
37 9911 2.3
38 8915 2.8
39 9452 8.8
40 9112 10.7
41 8472 13.9
42 8230 19.3
43 8384 19.5
44 8625 20.4
45 8221 15.3
46 8649 7.9
47 8625 8.3
48 10443 4.5
49 10357 3.2
50 8586 5.0
51 8892 6.6
52 8329 11.1
53 8101 12.8
54 7922 16.3
55 8120 17.4
56 7838 18.9
57 7735 15.8
58 8406 11.7
59 8209 6.4
60 9451 2.9
61 10041 4.7
62 9411 2.4
63 10405 7.2
64 8467 10.7
65 8464 13.4
66 8102 18.3
67 7627 18.4
68 7513 16.8
69 7510 16.6
70 8291 14.1
71 8064 6.1
72 9383 3.5
73 9706 1.7
74 8579 2.3
75 9474 4.5
76 8318 9.3
77 8213 14.2
78 8059 17.3
79 9111 23.0
80 7708 16.3
81 7680 18.4
82 8014 14.2
83 8007 9.1
84 8718 5.9
85 9486 7.2
86 9113 6.8
87 9025 8.0
88 8476 14.3
89 7952 14.6
90 7759 17.5
91 7835 17.2
92 7600 17.2
93 7651 14.1
94 8319 10.4
95 8812 6.8
96 8630 4.1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Temperatuur
9702.04 -96.54
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1049.13 -339.24 -52.18 177.66 2692.13
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9702.04 136.45 71.105 < 2e-16 ***
Temperatuur -96.54 11.00 -8.779 7.23e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 597 on 94 degrees of freedom
Multiple R-squared: 0.4505, Adjusted R-squared: 0.4447
F-statistic: 77.07 on 1 and 94 DF, p-value: 7.232e-14
> 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.9969123 0.006175379 0.003087689
[2,] 0.9957689 0.008462226 0.004231113
[3,] 0.9906417 0.018716647 0.009358323
[4,] 0.9881921 0.023615752 0.011807876
[5,] 0.9810185 0.037962986 0.018981493
[6,] 0.9747257 0.050548503 0.025274251
[7,] 0.9762765 0.047446965 0.023723482
[8,] 0.9772666 0.045466756 0.022733378
[9,] 0.9686620 0.062675978 0.031337989
[10,] 0.9855516 0.028896811 0.014448405
[11,] 0.9770748 0.045850382 0.022925191
[12,] 0.9696351 0.060729844 0.030364922
[13,] 0.9703073 0.059385360 0.029692680
[14,] 0.9560160 0.087967989 0.043983994
[15,] 0.9458558 0.108288315 0.054144158
[16,] 0.9332408 0.133518452 0.066759226
[17,] 0.9434984 0.113003294 0.056501647
[18,] 0.9287690 0.142461944 0.071230972
[19,] 0.9488184 0.102363271 0.051181636
[20,] 0.9286691 0.142661877 0.071330939
[21,] 0.9949930 0.010014091 0.005007045
[22,] 0.9923630 0.015274087 0.007637043
[23,] 0.9887195 0.022560998 0.011280499
[24,] 0.9832102 0.033579680 0.016789840
[25,] 0.9754876 0.049024824 0.024512412
[26,] 0.9649807 0.070038587 0.035019293
[27,] 0.9585406 0.082918867 0.041459433
[28,] 0.9494469 0.101106263 0.050553132
[29,] 0.9434890 0.113021987 0.056510993
[30,] 0.9248359 0.150328166 0.075164083
[31,] 0.9242693 0.151461359 0.075730680
[32,] 0.9103319 0.179336214 0.089668107
[33,] 0.8963075 0.207385036 0.103692518
[34,] 0.8938165 0.212367054 0.106183527
[35,] 0.8924627 0.215074661 0.107537331
[36,] 0.8788977 0.242204612 0.121102306
[37,] 0.8474902 0.305019695 0.152509848
[38,] 0.8262870 0.347425970 0.173712985
[39,] 0.8217126 0.356574723 0.178287361
[40,] 0.8666433 0.266713488 0.133356744
[41,] 0.8337107 0.332578565 0.166289282
[42,] 0.8054392 0.389121530 0.194560765
[43,] 0.7728584 0.454283281 0.227141640
[44,] 0.8874232 0.225153693 0.112576846
[45,] 0.9375728 0.124854376 0.062427188
[46,] 0.9391844 0.121631245 0.060815622
[47,] 0.9214625 0.157075073 0.078537537
[48,] 0.9033684 0.193263189 0.096631595
[49,] 0.8853715 0.229257075 0.114628537
[50,] 0.8572243 0.285551480 0.142775740
[51,] 0.8231541 0.353691767 0.176845883
[52,] 0.7817425 0.436514924 0.218257462
[53,] 0.7576764 0.484647164 0.242323582
[54,] 0.7109507 0.578098658 0.289049329
[55,] 0.7575427 0.484914693 0.242457346
[56,] 0.7109033 0.578193495 0.289096748
[57,] 0.7805010 0.438998034 0.219499017
[58,] 0.7376746 0.524650785 0.262325392
[59,] 0.9564945 0.087010975 0.043505488
[60,] 0.9405109 0.118978186 0.059489093
[61,] 0.9216460 0.156707958 0.078353979
[62,] 0.8993743 0.201251353 0.100625677
[63,] 0.8736970 0.252605961 0.126302980
[64,] 0.8677078 0.264584452 0.132292226
[65,] 0.8665435 0.266912973 0.133456487
[66,] 0.8260329 0.347934293 0.173967146
[67,] 0.8777455 0.244509078 0.122254539
[68,] 0.8536334 0.292733220 0.146366610
[69,] 0.8595762 0.280847569 0.140423784
[70,] 0.8597035 0.280592908 0.140296454
[71,] 0.8647266 0.270546860 0.135273430
[72,] 0.8297774 0.340445208 0.170222604
[73,] 0.7753393 0.449321442 0.224660721
[74,] 0.7104907 0.579018626 0.289509313
[75,] 0.9950431 0.009913890 0.004956945
[76,] 0.9911925 0.017614971 0.008807485
[77,] 0.9837888 0.032422438 0.016211219
[78,] 0.9711841 0.057631752 0.028815876
[79,] 0.9807335 0.038533025 0.019266513
[80,] 0.9708058 0.058388317 0.029194158
[81,] 0.9879514 0.024097140 0.012048570
[82,] 0.9847437 0.030512511 0.015256256
[83,] 0.9887080 0.022584005 0.011292003
[84,] 0.9962993 0.007401460 0.003700730
[85,] 0.9878964 0.024207274 0.012103637
[86,] 0.9662517 0.067496627 0.033748314
[87,] 0.9339993 0.132001305 0.066000653
> postscript(file="/var/www/rcomp/tmp/1a2kb1323605723.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/2rxg61323605723.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/3fc8y1323605723.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/4vl771323605723.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/5rf7s1323605723.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 = 96
Frequency = 1
1 2 3 4 5 6
2692.133812 36.565893 -228.582267 -271.335385 2.763337 116.886777
7 8 9 10 11 12
11.071977 327.355898 -414.693463 -245.100825 -313.347706 -298.397067
13 14 15 16 17 18
728.207892 -775.631627 39.491813 -318.421786 627.763337 114.689257
19 20 21 22 23 24
404.318858 429.590458 -720.631704 -338.816904 -835.471147 92.553572
25 26 27 28 29 30
1842.059733 -54.582267 223.997974 -8.261305 69.948536 26.849737
31 32 33 34 35 36
427.429978 384.664538 -472.853944 -72.335385 -559.767466 414.751092
37 38 39 40 41 42
431.010372 -516.718028 599.541174 442.973255 111.911496 391.244778
43 44 45 46 47 48
564.553418 892.442299 -3.928023 -290.347706 -275.730426 1175.405412
49 50 51 52 53 54
963.899252 -633.322987 -172.853867 -301.409465 -365.286024 -206.384823
55 56 57 58 59 60
97.812697 -39.372502 -441.656423 -166.483545 -875.162507 28.936292
61 62 63 64 65 66
792.714053 -59.335308 1398.072054 -202.026745 55.639896 166.701578
67 68 69 70 71 72
-298.644102 -567.113223 -589.421863 -49.779864 -1049.125467 18.862212
73 74 75 76 77 78
168.084451 -900.989628 206.405412 -486.187226 -118.125544 27.158377
79 80 81 82 83 84
1629.454620 -420.384823 -245.644102 -317.125544 -816.495866 -414.434107
85 86 87 88 89 90
479.072054 67.454773 95.306614 154.528776 -340.508264 -253.532983
91 92 93 94 95 96
-206.495943 -441.495943 -689.779864 -378.989705 -233.545227 -676.211868
> postscript(file="/var/www/rcomp/tmp/6ga1h1323605723.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 = 96
Frequency = 1
lag(myerror, k = 1) myerror
0 2692.133812 NA
1 36.565893 2692.133812
2 -228.582267 36.565893
3 -271.335385 -228.582267
4 2.763337 -271.335385
5 116.886777 2.763337
6 11.071977 116.886777
7 327.355898 11.071977
8 -414.693463 327.355898
9 -245.100825 -414.693463
10 -313.347706 -245.100825
11 -298.397067 -313.347706
12 728.207892 -298.397067
13 -775.631627 728.207892
14 39.491813 -775.631627
15 -318.421786 39.491813
16 627.763337 -318.421786
17 114.689257 627.763337
18 404.318858 114.689257
19 429.590458 404.318858
20 -720.631704 429.590458
21 -338.816904 -720.631704
22 -835.471147 -338.816904
23 92.553572 -835.471147
24 1842.059733 92.553572
25 -54.582267 1842.059733
26 223.997974 -54.582267
27 -8.261305 223.997974
28 69.948536 -8.261305
29 26.849737 69.948536
30 427.429978 26.849737
31 384.664538 427.429978
32 -472.853944 384.664538
33 -72.335385 -472.853944
34 -559.767466 -72.335385
35 414.751092 -559.767466
36 431.010372 414.751092
37 -516.718028 431.010372
38 599.541174 -516.718028
39 442.973255 599.541174
40 111.911496 442.973255
41 391.244778 111.911496
42 564.553418 391.244778
43 892.442299 564.553418
44 -3.928023 892.442299
45 -290.347706 -3.928023
46 -275.730426 -290.347706
47 1175.405412 -275.730426
48 963.899252 1175.405412
49 -633.322987 963.899252
50 -172.853867 -633.322987
51 -301.409465 -172.853867
52 -365.286024 -301.409465
53 -206.384823 -365.286024
54 97.812697 -206.384823
55 -39.372502 97.812697
56 -441.656423 -39.372502
57 -166.483545 -441.656423
58 -875.162507 -166.483545
59 28.936292 -875.162507
60 792.714053 28.936292
61 -59.335308 792.714053
62 1398.072054 -59.335308
63 -202.026745 1398.072054
64 55.639896 -202.026745
65 166.701578 55.639896
66 -298.644102 166.701578
67 -567.113223 -298.644102
68 -589.421863 -567.113223
69 -49.779864 -589.421863
70 -1049.125467 -49.779864
71 18.862212 -1049.125467
72 168.084451 18.862212
73 -900.989628 168.084451
74 206.405412 -900.989628
75 -486.187226 206.405412
76 -118.125544 -486.187226
77 27.158377 -118.125544
78 1629.454620 27.158377
79 -420.384823 1629.454620
80 -245.644102 -420.384823
81 -317.125544 -245.644102
82 -816.495866 -317.125544
83 -414.434107 -816.495866
84 479.072054 -414.434107
85 67.454773 479.072054
86 95.306614 67.454773
87 154.528776 95.306614
88 -340.508264 154.528776
89 -253.532983 -340.508264
90 -206.495943 -253.532983
91 -441.495943 -206.495943
92 -689.779864 -441.495943
93 -378.989705 -689.779864
94 -233.545227 -378.989705
95 -676.211868 -233.545227
96 NA -676.211868
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 36.565893 2692.133812
[2,] -228.582267 36.565893
[3,] -271.335385 -228.582267
[4,] 2.763337 -271.335385
[5,] 116.886777 2.763337
[6,] 11.071977 116.886777
[7,] 327.355898 11.071977
[8,] -414.693463 327.355898
[9,] -245.100825 -414.693463
[10,] -313.347706 -245.100825
[11,] -298.397067 -313.347706
[12,] 728.207892 -298.397067
[13,] -775.631627 728.207892
[14,] 39.491813 -775.631627
[15,] -318.421786 39.491813
[16,] 627.763337 -318.421786
[17,] 114.689257 627.763337
[18,] 404.318858 114.689257
[19,] 429.590458 404.318858
[20,] -720.631704 429.590458
[21,] -338.816904 -720.631704
[22,] -835.471147 -338.816904
[23,] 92.553572 -835.471147
[24,] 1842.059733 92.553572
[25,] -54.582267 1842.059733
[26,] 223.997974 -54.582267
[27,] -8.261305 223.997974
[28,] 69.948536 -8.261305
[29,] 26.849737 69.948536
[30,] 427.429978 26.849737
[31,] 384.664538 427.429978
[32,] -472.853944 384.664538
[33,] -72.335385 -472.853944
[34,] -559.767466 -72.335385
[35,] 414.751092 -559.767466
[36,] 431.010372 414.751092
[37,] -516.718028 431.010372
[38,] 599.541174 -516.718028
[39,] 442.973255 599.541174
[40,] 111.911496 442.973255
[41,] 391.244778 111.911496
[42,] 564.553418 391.244778
[43,] 892.442299 564.553418
[44,] -3.928023 892.442299
[45,] -290.347706 -3.928023
[46,] -275.730426 -290.347706
[47,] 1175.405412 -275.730426
[48,] 963.899252 1175.405412
[49,] -633.322987 963.899252
[50,] -172.853867 -633.322987
[51,] -301.409465 -172.853867
[52,] -365.286024 -301.409465
[53,] -206.384823 -365.286024
[54,] 97.812697 -206.384823
[55,] -39.372502 97.812697
[56,] -441.656423 -39.372502
[57,] -166.483545 -441.656423
[58,] -875.162507 -166.483545
[59,] 28.936292 -875.162507
[60,] 792.714053 28.936292
[61,] -59.335308 792.714053
[62,] 1398.072054 -59.335308
[63,] -202.026745 1398.072054
[64,] 55.639896 -202.026745
[65,] 166.701578 55.639896
[66,] -298.644102 166.701578
[67,] -567.113223 -298.644102
[68,] -589.421863 -567.113223
[69,] -49.779864 -589.421863
[70,] -1049.125467 -49.779864
[71,] 18.862212 -1049.125467
[72,] 168.084451 18.862212
[73,] -900.989628 168.084451
[74,] 206.405412 -900.989628
[75,] -486.187226 206.405412
[76,] -118.125544 -486.187226
[77,] 27.158377 -118.125544
[78,] 1629.454620 27.158377
[79,] -420.384823 1629.454620
[80,] -245.644102 -420.384823
[81,] -317.125544 -245.644102
[82,] -816.495866 -317.125544
[83,] -414.434107 -816.495866
[84,] 479.072054 -414.434107
[85,] 67.454773 479.072054
[86,] 95.306614 67.454773
[87,] 154.528776 95.306614
[88,] -340.508264 154.528776
[89,] -253.532983 -340.508264
[90,] -206.495943 -253.532983
[91,] -441.495943 -206.495943
[92,] -689.779864 -441.495943
[93,] -378.989705 -689.779864
[94,] -233.545227 -378.989705
[95,] -676.211868 -233.545227
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 36.565893 2692.133812
2 -228.582267 36.565893
3 -271.335385 -228.582267
4 2.763337 -271.335385
5 116.886777 2.763337
6 11.071977 116.886777
7 327.355898 11.071977
8 -414.693463 327.355898
9 -245.100825 -414.693463
10 -313.347706 -245.100825
11 -298.397067 -313.347706
12 728.207892 -298.397067
13 -775.631627 728.207892
14 39.491813 -775.631627
15 -318.421786 39.491813
16 627.763337 -318.421786
17 114.689257 627.763337
18 404.318858 114.689257
19 429.590458 404.318858
20 -720.631704 429.590458
21 -338.816904 -720.631704
22 -835.471147 -338.816904
23 92.553572 -835.471147
24 1842.059733 92.553572
25 -54.582267 1842.059733
26 223.997974 -54.582267
27 -8.261305 223.997974
28 69.948536 -8.261305
29 26.849737 69.948536
30 427.429978 26.849737
31 384.664538 427.429978
32 -472.853944 384.664538
33 -72.335385 -472.853944
34 -559.767466 -72.335385
35 414.751092 -559.767466
36 431.010372 414.751092
37 -516.718028 431.010372
38 599.541174 -516.718028
39 442.973255 599.541174
40 111.911496 442.973255
41 391.244778 111.911496
42 564.553418 391.244778
43 892.442299 564.553418
44 -3.928023 892.442299
45 -290.347706 -3.928023
46 -275.730426 -290.347706
47 1175.405412 -275.730426
48 963.899252 1175.405412
49 -633.322987 963.899252
50 -172.853867 -633.322987
51 -301.409465 -172.853867
52 -365.286024 -301.409465
53 -206.384823 -365.286024
54 97.812697 -206.384823
55 -39.372502 97.812697
56 -441.656423 -39.372502
57 -166.483545 -441.656423
58 -875.162507 -166.483545
59 28.936292 -875.162507
60 792.714053 28.936292
61 -59.335308 792.714053
62 1398.072054 -59.335308
63 -202.026745 1398.072054
64 55.639896 -202.026745
65 166.701578 55.639896
66 -298.644102 166.701578
67 -567.113223 -298.644102
68 -589.421863 -567.113223
69 -49.779864 -589.421863
70 -1049.125467 -49.779864
71 18.862212 -1049.125467
72 168.084451 18.862212
73 -900.989628 168.084451
74 206.405412 -900.989628
75 -486.187226 206.405412
76 -118.125544 -486.187226
77 27.158377 -118.125544
78 1629.454620 27.158377
79 -420.384823 1629.454620
80 -245.644102 -420.384823
81 -317.125544 -245.644102
82 -816.495866 -317.125544
83 -414.434107 -816.495866
84 479.072054 -414.434107
85 67.454773 479.072054
86 95.306614 67.454773
87 154.528776 95.306614
88 -340.508264 154.528776
89 -253.532983 -340.508264
90 -206.495943 -253.532983
91 -441.495943 -206.495943
92 -689.779864 -441.495943
93 -378.989705 -689.779864
94 -233.545227 -378.989705
95 -676.211868 -233.545227
> 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/7j8rh1323605723.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/8oqry1323605723.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/9tc9l1323605723.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/109pso1323605723.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/1110nw1323605723.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/12lv4d1323605723.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/13ab5q1323605723.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/14l6851323605723.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/15e0m21323605723.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/16qe0k1323605723.tab")
+ }
>
> try(system("convert tmp/1a2kb1323605723.ps tmp/1a2kb1323605723.png",intern=TRUE))
character(0)
> try(system("convert tmp/2rxg61323605723.ps tmp/2rxg61323605723.png",intern=TRUE))
character(0)
> try(system("convert tmp/3fc8y1323605723.ps tmp/3fc8y1323605723.png",intern=TRUE))
character(0)
> try(system("convert tmp/4vl771323605723.ps tmp/4vl771323605723.png",intern=TRUE))
character(0)
> try(system("convert tmp/5rf7s1323605723.ps tmp/5rf7s1323605723.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ga1h1323605723.ps tmp/6ga1h1323605723.png",intern=TRUE))
character(0)
> try(system("convert tmp/7j8rh1323605723.ps tmp/7j8rh1323605723.png",intern=TRUE))
character(0)
> try(system("convert tmp/8oqry1323605723.ps tmp/8oqry1323605723.png",intern=TRUE))
character(0)
> try(system("convert tmp/9tc9l1323605723.ps tmp/9tc9l1323605723.png",intern=TRUE))
character(0)
> try(system("convert tmp/109pso1323605723.ps tmp/109pso1323605723.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.104 0.616 4.809