R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
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(2649.2
+ ,31077
+ ,0
+ ,2579.4
+ ,31293
+ ,0
+ ,2504.6
+ ,30236
+ ,0
+ ,2462.3
+ ,30160
+ ,0
+ ,2467.4
+ ,32436
+ ,0
+ ,2446.7
+ ,30695
+ ,0
+ ,2656.3
+ ,27525
+ ,0
+ ,2626.2
+ ,26434
+ ,0
+ ,2482.6
+ ,25739
+ ,0
+ ,2539.9
+ ,25204
+ ,0
+ ,2502.7
+ ,24977
+ ,0
+ ,2466.9
+ ,24320
+ ,0
+ ,2513.2
+ ,22680
+ ,1
+ ,2443.3
+ ,22052
+ ,1
+ ,2293.4
+ ,21467
+ ,1
+ ,2070.8
+ ,21383
+ ,1
+ ,2029.6
+ ,21777
+ ,1
+ ,2052
+ ,21928
+ ,1
+ ,1864.4
+ ,21814
+ ,1
+ ,1670.1
+ ,22937
+ ,1
+ ,1811
+ ,23595
+ ,1
+ ,1905.4
+ ,20830
+ ,1
+ ,1862.8
+ ,19650
+ ,1
+ ,2014.5
+ ,19195
+ ,1
+ ,2197.8
+ ,19644
+ ,1
+ ,2962.3
+ ,18483
+ ,0
+ ,3047
+ ,18079
+ ,0
+ ,3032.6
+ ,19178
+ ,0
+ ,3504.4
+ ,18391
+ ,0
+ ,3801.1
+ ,18441
+ ,0
+ ,3857.6
+ ,18584
+ ,0
+ ,3674.4
+ ,20108
+ ,0
+ ,3721
+ ,20148
+ ,0
+ ,3844.5
+ ,19394
+ ,0
+ ,4116.7
+ ,17745
+ ,0
+ ,4105.2
+ ,17696
+ ,0
+ ,4435.2
+ ,17032
+ ,0
+ ,4296.5
+ ,16438
+ ,0
+ ,4202.5
+ ,15683
+ ,0
+ ,4562.8
+ ,15594
+ ,0
+ ,4621.4
+ ,15713
+ ,0
+ ,4697
+ ,15937
+ ,0
+ ,4591.3
+ ,16171
+ ,0
+ ,4357
+ ,15928
+ ,0
+ ,4502.6
+ ,16348
+ ,0
+ ,4443.9
+ ,15579
+ ,0
+ ,4290.9
+ ,15305
+ ,0
+ ,4199.8
+ ,15648
+ ,0
+ ,4138.5
+ ,14954
+ ,0
+ ,3970.1
+ ,15137
+ ,0
+ ,3862.3
+ ,15839
+ ,0
+ ,3701.6
+ ,16050
+ ,0
+ ,3570.12
+ ,15168
+ ,0
+ ,3801.06
+ ,17064
+ ,0
+ ,3895.51
+ ,16005
+ ,0
+ ,3917.96
+ ,14886
+ ,0
+ ,3813.06
+ ,14931
+ ,0
+ ,3667.03
+ ,14544
+ ,0
+ ,3494.17
+ ,13812
+ ,0
+ ,3364
+ ,13031
+ ,0
+ ,3295.3
+ ,12574
+ ,0
+ ,3277.0
+ ,11964
+ ,0
+ ,3257.2
+ ,11451
+ ,0
+ ,3161.7
+ ,11346
+ ,0
+ ,3097.3
+ ,11353
+ ,0
+ ,3061.3
+ ,10702
+ ,0
+ ,3119.3
+ ,10646
+ ,0
+ ,3106.22
+ ,10556
+ ,0
+ ,3080.58
+ ,10463
+ ,0
+ ,2981.85
+ ,10407
+ ,0)
+ ,dim=c(3
+ ,70)
+ ,dimnames=list(c('Bel20'
+ ,'Goudprijs'
+ ,'Crisis')
+ ,1:70))
> y <- array(NA,dim=c(3,70),dimnames=list(c('Bel20','Goudprijs','Crisis'),1:70))
> 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
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Bel20 Goudprijs Crisis
1 2649.20 31077 0
2 2579.40 31293 0
3 2504.60 30236 0
4 2462.30 30160 0
5 2467.40 32436 0
6 2446.70 30695 0
7 2656.30 27525 0
8 2626.20 26434 0
9 2482.60 25739 0
10 2539.90 25204 0
11 2502.70 24977 0
12 2466.90 24320 0
13 2513.20 22680 1
14 2443.30 22052 1
15 2293.40 21467 1
16 2070.80 21383 1
17 2029.60 21777 1
18 2052.00 21928 1
19 1864.40 21814 1
20 1670.10 22937 1
21 1811.00 23595 1
22 1905.40 20830 1
23 1862.80 19650 1
24 2014.50 19195 1
25 2197.80 19644 1
26 2962.30 18483 0
27 3047.00 18079 0
28 3032.60 19178 0
29 3504.40 18391 0
30 3801.10 18441 0
31 3857.60 18584 0
32 3674.40 20108 0
33 3721.00 20148 0
34 3844.50 19394 0
35 4116.70 17745 0
36 4105.20 17696 0
37 4435.20 17032 0
38 4296.50 16438 0
39 4202.50 15683 0
40 4562.80 15594 0
41 4621.40 15713 0
42 4697.00 15937 0
43 4591.30 16171 0
44 4357.00 15928 0
45 4502.60 16348 0
46 4443.90 15579 0
47 4290.90 15305 0
48 4199.80 15648 0
49 4138.50 14954 0
50 3970.10 15137 0
51 3862.30 15839 0
52 3701.60 16050 0
53 3570.12 15168 0
54 3801.06 17064 0
55 3895.51 16005 0
56 3917.96 14886 0
57 3813.06 14931 0
58 3667.03 14544 0
59 3494.17 13812 0
60 3364.00 13031 0
61 3295.30 12574 0
62 3277.00 11964 0
63 3257.20 11451 0
64 3161.70 11346 0
65 3097.30 11353 0
66 3061.30 10702 0
67 3119.30 10646 0
68 3106.22 10556 0
69 3080.58 10463 0
70 2981.85 10407 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Goudprijs Crisis
4613.2218 -0.0612 -1244.0672
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-994.51 -397.60 -31.53 408.99 1059.06
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.613e+03 2.283e+02 20.205 < 2e-16 ***
Goudprijs -6.120e-02 1.199e-02 -5.106 2.93e-06 ***
Crisis -1.244e+03 1.717e+02 -7.247 5.47e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 542.8 on 67 degrees of freedom
Multiple R-squared: 0.6028, Adjusted R-squared: 0.591
F-statistic: 50.84 on 2 and 67 DF, p-value: 3.684e-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,] 8.343082e-03 1.668616e-02 9.916569e-01
[2,] 1.818705e-03 3.637410e-03 9.981813e-01
[3,] 2.870938e-04 5.741877e-04 9.997129e-01
[4,] 1.308132e-04 2.616263e-04 9.998692e-01
[5,] 2.831672e-05 5.663345e-05 9.999717e-01
[6,] 8.960722e-06 1.792144e-05 9.999910e-01
[7,] 6.170766e-06 1.234153e-05 9.999938e-01
[8,] 1.141329e-06 2.282657e-06 9.999989e-01
[9,] 2.477402e-07 4.954805e-07 9.999998e-01
[10,] 2.255083e-07 4.510167e-07 9.999998e-01
[11,] 2.798179e-06 5.596359e-06 9.999972e-01
[12,] 5.966764e-06 1.193353e-05 9.999940e-01
[13,] 4.415200e-06 8.830400e-06 9.999956e-01
[14,] 1.370926e-05 2.741851e-05 9.999863e-01
[15,] 1.214120e-04 2.428239e-04 9.998786e-01
[16,] 1.252004e-04 2.504008e-04 9.998748e-01
[17,] 6.998534e-05 1.399707e-04 9.999300e-01
[18,] 4.233559e-05 8.467119e-05 9.999577e-01
[19,] 1.673335e-05 3.346669e-05 9.999833e-01
[20,] 8.031001e-06 1.606200e-05 9.999920e-01
[21,] 2.424689e-05 4.849377e-05 9.999758e-01
[22,] 5.179204e-05 1.035841e-04 9.999482e-01
[23,] 1.892060e-04 3.784121e-04 9.998108e-01
[24,] 1.612037e-03 3.224073e-03 9.983880e-01
[25,] 1.354144e-02 2.708287e-02 9.864586e-01
[26,] 4.220069e-02 8.440137e-02 9.577993e-01
[27,] 1.221215e-01 2.442430e-01 8.778785e-01
[28,] 3.584274e-01 7.168548e-01 6.415726e-01
[29,] 6.852266e-01 6.295469e-01 3.147734e-01
[30,] 7.873787e-01 4.252425e-01 2.126213e-01
[31,] 8.558057e-01 2.883886e-01 1.441943e-01
[32,] 8.949470e-01 2.101059e-01 1.050530e-01
[33,] 8.870476e-01 2.259048e-01 1.129524e-01
[34,] 8.588471e-01 2.823059e-01 1.411529e-01
[35,] 9.083411e-01 1.833179e-01 9.165895e-02
[36,] 9.534032e-01 9.319365e-02 4.659683e-02
[37,] 9.856870e-01 2.862593e-02 1.431296e-02
[38,] 9.935630e-01 1.287405e-02 6.437026e-03
[39,] 9.936912e-01 1.261768e-02 6.308838e-03
[40,] 9.965920e-01 6.816096e-03 3.408048e-03
[41,] 9.993214e-01 1.357147e-03 6.785734e-04
[42,] 9.998704e-01 2.592372e-04 1.296186e-04
[43,] 9.999534e-01 9.310237e-05 4.655118e-05
[44,] 9.999982e-01 3.661097e-06 1.830548e-06
[45,] 9.999995e-01 9.917948e-07 4.958974e-07
[46,] 9.999987e-01 2.653119e-06 1.326559e-06
[47,] 9.999979e-01 4.233684e-06 2.116842e-06
[48,] 9.999981e-01 3.823993e-06 1.911996e-06
[49,] 9.999998e-01 3.541120e-07 1.770560e-07
[50,] 9.999993e-01 1.471433e-06 7.357164e-07
[51,] 9.999998e-01 3.162103e-07 1.581052e-07
[52,] 9.999999e-01 2.323608e-07 1.161804e-07
[53,] 9.999997e-01 5.610374e-07 2.805187e-07
[54,] 9.999985e-01 3.007153e-06 1.503576e-06
[55,] 9.999922e-01 1.559608e-05 7.798039e-06
[56,] 9.999663e-01 6.746387e-05 3.373194e-05
[57,] 9.998065e-01 3.870964e-04 1.935482e-04
[58,] 9.995745e-01 8.509873e-04 4.254936e-04
[59,] 9.972368e-01 5.526499e-03 2.763249e-03
> postscript(file="/var/www/html/rcomp/tmp/1kq6z1290542319.ps",horizontal=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/html/rcomp/tmp/2kq6z1290542319.ps",horizontal=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/html/rcomp/tmp/3vzn21290542319.ps",horizontal=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/html/rcomp/tmp/4vzn21290542319.ps",horizontal=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/html/rcomp/tmp/5vzn21290542319.ps",horizontal=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 = 70
Frequency = 1
1 2 3 4 5 6
-62.237868 -118.819561 -258.303591 -305.254477 -160.872687 -288.114689
7 8 9 10 11 12
-272.505582 -369.370271 -555.501397 -530.941185 -582.032647 -658.038330
13 14 15 16 17 18
531.967685 423.636682 237.937101 10.196648 -6.892181 24.748395
19 20 21 22 23 24
-169.827934 -295.404977 -114.238098 -189.044666 -303.855787 -179.999905
25 26 27 28 29 30
30.777038 -519.838597 -459.861726 -407.007470 16.631384 316.391177
31 32 33 34 35 36
381.642186 291.704685 340.752519 418.110837 589.398855 574.900257
37 38 39 40 41 42
864.266203 689.215859 549.012980 903.866548 969.748856 1059.056730
43 44 45 46 47 48
967.676563 718.505967 889.808231 784.048610 614.280943 544.171125
49 50 51 52 53 54
440.401195 283.200038 218.359535 70.571863 -114.882890 232.084471
55 56 57 58 59 60
261.728049 215.699876 113.553690 -56.159110 -273.814484 -451.778454
61 62 63 64 65 66
-548.444965 -604.074443 -655.267922 -757.193488 -821.165117 -897.003625
67 68 69 70
-842.430593 -861.018221 -892.349437 -994.506405
> postscript(file="/var/www/html/rcomp/tmp/66q451290542319.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 70
Frequency = 1
lag(myerror, k = 1) myerror
0 -62.237868 NA
1 -118.819561 -62.237868
2 -258.303591 -118.819561
3 -305.254477 -258.303591
4 -160.872687 -305.254477
5 -288.114689 -160.872687
6 -272.505582 -288.114689
7 -369.370271 -272.505582
8 -555.501397 -369.370271
9 -530.941185 -555.501397
10 -582.032647 -530.941185
11 -658.038330 -582.032647
12 531.967685 -658.038330
13 423.636682 531.967685
14 237.937101 423.636682
15 10.196648 237.937101
16 -6.892181 10.196648
17 24.748395 -6.892181
18 -169.827934 24.748395
19 -295.404977 -169.827934
20 -114.238098 -295.404977
21 -189.044666 -114.238098
22 -303.855787 -189.044666
23 -179.999905 -303.855787
24 30.777038 -179.999905
25 -519.838597 30.777038
26 -459.861726 -519.838597
27 -407.007470 -459.861726
28 16.631384 -407.007470
29 316.391177 16.631384
30 381.642186 316.391177
31 291.704685 381.642186
32 340.752519 291.704685
33 418.110837 340.752519
34 589.398855 418.110837
35 574.900257 589.398855
36 864.266203 574.900257
37 689.215859 864.266203
38 549.012980 689.215859
39 903.866548 549.012980
40 969.748856 903.866548
41 1059.056730 969.748856
42 967.676563 1059.056730
43 718.505967 967.676563
44 889.808231 718.505967
45 784.048610 889.808231
46 614.280943 784.048610
47 544.171125 614.280943
48 440.401195 544.171125
49 283.200038 440.401195
50 218.359535 283.200038
51 70.571863 218.359535
52 -114.882890 70.571863
53 232.084471 -114.882890
54 261.728049 232.084471
55 215.699876 261.728049
56 113.553690 215.699876
57 -56.159110 113.553690
58 -273.814484 -56.159110
59 -451.778454 -273.814484
60 -548.444965 -451.778454
61 -604.074443 -548.444965
62 -655.267922 -604.074443
63 -757.193488 -655.267922
64 -821.165117 -757.193488
65 -897.003625 -821.165117
66 -842.430593 -897.003625
67 -861.018221 -842.430593
68 -892.349437 -861.018221
69 -994.506405 -892.349437
70 NA -994.506405
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -118.819561 -62.237868
[2,] -258.303591 -118.819561
[3,] -305.254477 -258.303591
[4,] -160.872687 -305.254477
[5,] -288.114689 -160.872687
[6,] -272.505582 -288.114689
[7,] -369.370271 -272.505582
[8,] -555.501397 -369.370271
[9,] -530.941185 -555.501397
[10,] -582.032647 -530.941185
[11,] -658.038330 -582.032647
[12,] 531.967685 -658.038330
[13,] 423.636682 531.967685
[14,] 237.937101 423.636682
[15,] 10.196648 237.937101
[16,] -6.892181 10.196648
[17,] 24.748395 -6.892181
[18,] -169.827934 24.748395
[19,] -295.404977 -169.827934
[20,] -114.238098 -295.404977
[21,] -189.044666 -114.238098
[22,] -303.855787 -189.044666
[23,] -179.999905 -303.855787
[24,] 30.777038 -179.999905
[25,] -519.838597 30.777038
[26,] -459.861726 -519.838597
[27,] -407.007470 -459.861726
[28,] 16.631384 -407.007470
[29,] 316.391177 16.631384
[30,] 381.642186 316.391177
[31,] 291.704685 381.642186
[32,] 340.752519 291.704685
[33,] 418.110837 340.752519
[34,] 589.398855 418.110837
[35,] 574.900257 589.398855
[36,] 864.266203 574.900257
[37,] 689.215859 864.266203
[38,] 549.012980 689.215859
[39,] 903.866548 549.012980
[40,] 969.748856 903.866548
[41,] 1059.056730 969.748856
[42,] 967.676563 1059.056730
[43,] 718.505967 967.676563
[44,] 889.808231 718.505967
[45,] 784.048610 889.808231
[46,] 614.280943 784.048610
[47,] 544.171125 614.280943
[48,] 440.401195 544.171125
[49,] 283.200038 440.401195
[50,] 218.359535 283.200038
[51,] 70.571863 218.359535
[52,] -114.882890 70.571863
[53,] 232.084471 -114.882890
[54,] 261.728049 232.084471
[55,] 215.699876 261.728049
[56,] 113.553690 215.699876
[57,] -56.159110 113.553690
[58,] -273.814484 -56.159110
[59,] -451.778454 -273.814484
[60,] -548.444965 -451.778454
[61,] -604.074443 -548.444965
[62,] -655.267922 -604.074443
[63,] -757.193488 -655.267922
[64,] -821.165117 -757.193488
[65,] -897.003625 -821.165117
[66,] -842.430593 -897.003625
[67,] -861.018221 -842.430593
[68,] -892.349437 -861.018221
[69,] -994.506405 -892.349437
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -118.819561 -62.237868
2 -258.303591 -118.819561
3 -305.254477 -258.303591
4 -160.872687 -305.254477
5 -288.114689 -160.872687
6 -272.505582 -288.114689
7 -369.370271 -272.505582
8 -555.501397 -369.370271
9 -530.941185 -555.501397
10 -582.032647 -530.941185
11 -658.038330 -582.032647
12 531.967685 -658.038330
13 423.636682 531.967685
14 237.937101 423.636682
15 10.196648 237.937101
16 -6.892181 10.196648
17 24.748395 -6.892181
18 -169.827934 24.748395
19 -295.404977 -169.827934
20 -114.238098 -295.404977
21 -189.044666 -114.238098
22 -303.855787 -189.044666
23 -179.999905 -303.855787
24 30.777038 -179.999905
25 -519.838597 30.777038
26 -459.861726 -519.838597
27 -407.007470 -459.861726
28 16.631384 -407.007470
29 316.391177 16.631384
30 381.642186 316.391177
31 291.704685 381.642186
32 340.752519 291.704685
33 418.110837 340.752519
34 589.398855 418.110837
35 574.900257 589.398855
36 864.266203 574.900257
37 689.215859 864.266203
38 549.012980 689.215859
39 903.866548 549.012980
40 969.748856 903.866548
41 1059.056730 969.748856
42 967.676563 1059.056730
43 718.505967 967.676563
44 889.808231 718.505967
45 784.048610 889.808231
46 614.280943 784.048610
47 544.171125 614.280943
48 440.401195 544.171125
49 283.200038 440.401195
50 218.359535 283.200038
51 70.571863 218.359535
52 -114.882890 70.571863
53 232.084471 -114.882890
54 261.728049 232.084471
55 215.699876 261.728049
56 113.553690 215.699876
57 -56.159110 113.553690
58 -273.814484 -56.159110
59 -451.778454 -273.814484
60 -548.444965 -451.778454
61 -604.074443 -548.444965
62 -655.267922 -604.074443
63 -757.193488 -655.267922
64 -821.165117 -757.193488
65 -897.003625 -821.165117
66 -842.430593 -897.003625
67 -861.018221 -842.430593
68 -892.349437 -861.018221
69 -994.506405 -892.349437
> 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/html/rcomp/tmp/7y0lq1290542319.ps",horizontal=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/html/rcomp/tmp/8y0lq1290542319.ps",horizontal=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/html/rcomp/tmp/9y0lq1290542319.ps",horizontal=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/html/rcomp/tmp/10rrlt1290542319.ps",horizontal=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/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/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/html/rcomp/tmp/11ds1h1290542319.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/html/rcomp/tmp/12ys041290542319.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/html/rcomp/tmp/13nbfg1290542319.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/html/rcomp/tmp/14xkwj1290542319.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/html/rcomp/tmp/151lvp1290542319.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/html/rcomp/tmp/16xuay1290542319.tab")
+ }
>
> try(system("convert tmp/1kq6z1290542319.ps tmp/1kq6z1290542319.png",intern=TRUE))
character(0)
> try(system("convert tmp/2kq6z1290542319.ps tmp/2kq6z1290542319.png",intern=TRUE))
character(0)
> try(system("convert tmp/3vzn21290542319.ps tmp/3vzn21290542319.png",intern=TRUE))
character(0)
> try(system("convert tmp/4vzn21290542319.ps tmp/4vzn21290542319.png",intern=TRUE))
character(0)
> try(system("convert tmp/5vzn21290542319.ps tmp/5vzn21290542319.png",intern=TRUE))
character(0)
> try(system("convert tmp/66q451290542319.ps tmp/66q451290542319.png",intern=TRUE))
character(0)
> try(system("convert tmp/7y0lq1290542319.ps tmp/7y0lq1290542319.png",intern=TRUE))
character(0)
> try(system("convert tmp/8y0lq1290542319.ps tmp/8y0lq1290542319.png",intern=TRUE))
character(0)
> try(system("convert tmp/9y0lq1290542319.ps tmp/9y0lq1290542319.png",intern=TRUE))
character(0)
> try(system("convert tmp/10rrlt1290542319.ps tmp/10rrlt1290542319.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.505 1.568 5.859