R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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.
Natural language support but running in an English locale
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 = '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 t
1 2649.20 31077 0 1
2 2579.40 31293 0 2
3 2504.60 30236 0 3
4 2462.30 30160 0 4
5 2467.40 32436 0 5
6 2446.70 30695 0 6
7 2656.30 27525 0 7
8 2626.20 26434 0 8
9 2482.60 25739 0 9
10 2539.90 25204 0 10
11 2502.70 24977 0 11
12 2466.90 24320 0 12
13 2513.20 22680 1 13
14 2443.30 22052 1 14
15 2293.40 21467 1 15
16 2070.80 21383 1 16
17 2029.60 21777 1 17
18 2052.00 21928 1 18
19 1864.40 21814 1 19
20 1670.10 22937 1 20
21 1811.00 23595 1 21
22 1905.40 20830 1 22
23 1862.80 19650 1 23
24 2014.50 19195 1 24
25 2197.80 19644 1 25
26 2962.30 18483 0 26
27 3047.00 18079 0 27
28 3032.60 19178 0 28
29 3504.40 18391 0 29
30 3801.10 18441 0 30
31 3857.60 18584 0 31
32 3674.40 20108 0 32
33 3721.00 20148 0 33
34 3844.50 19394 0 34
35 4116.70 17745 0 35
36 4105.20 17696 0 36
37 4435.20 17032 0 37
38 4296.50 16438 0 38
39 4202.50 15683 0 39
40 4562.80 15594 0 40
41 4621.40 15713 0 41
42 4697.00 15937 0 42
43 4591.30 16171 0 43
44 4357.00 15928 0 44
45 4502.60 16348 0 45
46 4443.90 15579 0 46
47 4290.90 15305 0 47
48 4199.80 15648 0 48
49 4138.50 14954 0 49
50 3970.10 15137 0 50
51 3862.30 15839 0 51
52 3701.60 16050 0 52
53 3570.12 15168 0 53
54 3801.06 17064 0 54
55 3895.51 16005 0 55
56 3917.96 14886 0 56
57 3813.06 14931 0 57
58 3667.03 14544 0 58
59 3494.17 13812 0 59
60 3364.00 13031 0 60
61 3295.30 12574 0 61
62 3277.00 11964 0 62
63 3257.20 11451 0 63
64 3161.70 11346 0 64
65 3097.30 11353 0 65
66 3061.30 10702 0 66
67 3119.30 10646 0 67
68 3106.22 10556 0 68
69 3080.58 10463 0 69
70 2981.85 10407 0 70
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Goudprijs Crisis t
7157.5498 -0.1459 -1480.7553 -25.7971
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-851.55 -429.04 38.88 389.39 948.08
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.158e+03 1.150e+03 6.224 3.79e-08 ***
Goudprijs -1.459e-01 3.933e-02 -3.710 0.000428 ***
Crisis -1.481e+03 1.970e+02 -7.517 1.93e-10 ***
t -2.580e+01 1.144e+01 -2.255 0.027484 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 527 on 66 degrees of freedom
Multiple R-squared: 0.6312, Adjusted R-squared: 0.6145
F-statistic: 37.66 on 3 and 66 DF, p-value: 2.654e-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,] 5.800093e-03 1.160019e-02 9.941999e-01
[2,] 7.853531e-04 1.570706e-03 9.992146e-01
[3,] 2.828283e-04 5.656565e-04 9.997172e-01
[4,] 4.781839e-05 9.563678e-05 9.999522e-01
[5,] 9.957801e-06 1.991560e-05 9.999900e-01
[6,] 3.963871e-06 7.927742e-06 9.999960e-01
[7,] 5.859381e-07 1.171876e-06 9.999994e-01
[8,] 1.035375e-07 2.070749e-07 9.999999e-01
[9,] 1.092538e-07 2.185077e-07 9.999999e-01
[10,] 1.231271e-06 2.462541e-06 9.999988e-01
[11,] 8.130403e-07 1.626081e-06 9.999992e-01
[12,] 1.873481e-07 3.746962e-07 9.999998e-01
[13,] 8.973008e-08 1.794602e-07 9.999999e-01
[14,] 3.809700e-08 7.619399e-08 1.000000e+00
[15,] 2.087936e-08 4.175872e-08 1.000000e+00
[16,] 5.432558e-09 1.086512e-08 1.000000e+00
[17,] 1.197393e-09 2.394786e-09 1.000000e+00
[18,] 7.554235e-10 1.510847e-09 1.000000e+00
[19,] 1.677307e-08 3.354613e-08 1.000000e+00
[20,] 3.571394e-06 7.142789e-06 9.999964e-01
[21,] 2.897564e-05 5.795128e-05 9.999710e-01
[22,] 3.471383e-04 6.942766e-04 9.996529e-01
[23,] 8.920675e-03 1.784135e-02 9.910793e-01
[24,] 8.828805e-02 1.765761e-01 9.117120e-01
[25,] 2.490482e-01 4.980964e-01 7.509518e-01
[26,] 4.434805e-01 8.869611e-01 5.565195e-01
[27,] 7.131806e-01 5.736388e-01 2.868194e-01
[28,] 9.339399e-01 1.321202e-01 6.606009e-02
[29,] 9.849822e-01 3.003559e-02 1.501779e-02
[30,] 9.986141e-01 2.771773e-03 1.385887e-03
[31,] 9.994746e-01 1.050815e-03 5.254074e-04
[32,] 9.998477e-01 3.045430e-04 1.522715e-04
[33,] 9.999925e-01 1.500897e-05 7.504485e-06
[34,] 9.999924e-01 1.514050e-05 7.570251e-06
[35,] 9.999881e-01 2.375225e-05 1.187613e-05
[36,] 9.999867e-01 2.668547e-05 1.334274e-05
[37,] 9.999780e-01 4.391211e-05 2.195606e-05
[38,] 9.999492e-01 1.016421e-04 5.082106e-05
[39,] 9.999323e-01 1.354033e-04 6.770165e-05
[40,] 9.999406e-01 1.187157e-04 5.935784e-05
[41,] 9.999372e-01 1.256692e-04 6.283459e-05
[42,] 9.999365e-01 1.270641e-04 6.353207e-05
[43,] 9.999801e-01 3.973950e-05 1.986975e-05
[44,] 9.999917e-01 1.652214e-05 8.261071e-06
[45,] 9.999878e-01 2.444687e-05 1.222343e-05
[46,] 9.999878e-01 2.435244e-05 1.217622e-05
[47,] 9.999949e-01 1.024577e-05 5.122885e-06
[48,] 9.999991e-01 1.795903e-06 8.979516e-07
[49,] 9.999967e-01 6.604415e-06 3.302208e-06
[50,] 9.999992e-01 1.685214e-06 8.426069e-07
[51,] 9.999993e-01 1.462253e-06 7.311267e-07
[52,] 9.999989e-01 2.268219e-06 1.134109e-06
[53,] 9.999963e-01 7.456531e-06 3.728266e-06
[54,] 9.999793e-01 4.146254e-05 2.073127e-05
[55,] 9.998681e-01 2.638210e-04 1.319105e-04
[56,] 9.993165e-01 1.366968e-03 6.834842e-04
[57,] 9.976540e-01 4.691912e-03 2.345956e-03
> postscript(file="/var/www/html/freestat/rcomp/tmp/1aj2w1291108641.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/html/freestat/rcomp/tmp/2aj2w1291108641.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/html/freestat/rcomp/tmp/3aj2w1291108641.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/html/freestat/rcomp/tmp/4la1z1291108641.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/html/freestat/rcomp/tmp/5la1z1291108641.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 = 70
Frequency = 1
1 2 3 4 5 6
51.474372 38.985168 -164.230340 -191.821336 171.136347 -77.772400
7 8 9 10 11 12
-304.867312 -468.343303 -687.544262 -682.501773 -727.023147 -832.880037
13 14 15 16 17 18
480.702099 344.976209 135.523871 -73.534298 -31.453916 38.773604
19 20 21 22 23 24
-139.661461 -144.322497 118.374573 -164.832238 -353.793021 -242.678808
25 26 27 28 29 30
31.925883 -827.918197 -776.363260 -604.625813 -221.849254 107.942715
31 32 33 34 35 36
211.103063 276.046542 354.279546 393.570691 450.984427 458.132638
37 38 39 40 41 42
717.054472 517.489065 339.134314 712.246663 814.005494 948.083462
43 44 45 46 47 48
902.320396 658.364677 891.038368 745.941065 578.762553 563.502211
49 50 51 52 53 54
426.747149 310.843359 331.259876 227.141189 -7.222423 526.134572
55 56 57 58 59 60
491.877270 376.866176 304.328662 127.633840 -126.225291 -344.543352
61 62 63 64 65 66
-454.120932 -535.620684 -604.468471 -689.490467 -727.072049 -832.253560
67 68 69 70
-756.626624 -757.040172 -770.451409 -851.554474
> postscript(file="/var/www/html/freestat/rcomp/tmp/6la1z1291108641.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 = 70
Frequency = 1
lag(myerror, k = 1) myerror
0 51.474372 NA
1 38.985168 51.474372
2 -164.230340 38.985168
3 -191.821336 -164.230340
4 171.136347 -191.821336
5 -77.772400 171.136347
6 -304.867312 -77.772400
7 -468.343303 -304.867312
8 -687.544262 -468.343303
9 -682.501773 -687.544262
10 -727.023147 -682.501773
11 -832.880037 -727.023147
12 480.702099 -832.880037
13 344.976209 480.702099
14 135.523871 344.976209
15 -73.534298 135.523871
16 -31.453916 -73.534298
17 38.773604 -31.453916
18 -139.661461 38.773604
19 -144.322497 -139.661461
20 118.374573 -144.322497
21 -164.832238 118.374573
22 -353.793021 -164.832238
23 -242.678808 -353.793021
24 31.925883 -242.678808
25 -827.918197 31.925883
26 -776.363260 -827.918197
27 -604.625813 -776.363260
28 -221.849254 -604.625813
29 107.942715 -221.849254
30 211.103063 107.942715
31 276.046542 211.103063
32 354.279546 276.046542
33 393.570691 354.279546
34 450.984427 393.570691
35 458.132638 450.984427
36 717.054472 458.132638
37 517.489065 717.054472
38 339.134314 517.489065
39 712.246663 339.134314
40 814.005494 712.246663
41 948.083462 814.005494
42 902.320396 948.083462
43 658.364677 902.320396
44 891.038368 658.364677
45 745.941065 891.038368
46 578.762553 745.941065
47 563.502211 578.762553
48 426.747149 563.502211
49 310.843359 426.747149
50 331.259876 310.843359
51 227.141189 331.259876
52 -7.222423 227.141189
53 526.134572 -7.222423
54 491.877270 526.134572
55 376.866176 491.877270
56 304.328662 376.866176
57 127.633840 304.328662
58 -126.225291 127.633840
59 -344.543352 -126.225291
60 -454.120932 -344.543352
61 -535.620684 -454.120932
62 -604.468471 -535.620684
63 -689.490467 -604.468471
64 -727.072049 -689.490467
65 -832.253560 -727.072049
66 -756.626624 -832.253560
67 -757.040172 -756.626624
68 -770.451409 -757.040172
69 -851.554474 -770.451409
70 NA -851.554474
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 38.985168 51.474372
[2,] -164.230340 38.985168
[3,] -191.821336 -164.230340
[4,] 171.136347 -191.821336
[5,] -77.772400 171.136347
[6,] -304.867312 -77.772400
[7,] -468.343303 -304.867312
[8,] -687.544262 -468.343303
[9,] -682.501773 -687.544262
[10,] -727.023147 -682.501773
[11,] -832.880037 -727.023147
[12,] 480.702099 -832.880037
[13,] 344.976209 480.702099
[14,] 135.523871 344.976209
[15,] -73.534298 135.523871
[16,] -31.453916 -73.534298
[17,] 38.773604 -31.453916
[18,] -139.661461 38.773604
[19,] -144.322497 -139.661461
[20,] 118.374573 -144.322497
[21,] -164.832238 118.374573
[22,] -353.793021 -164.832238
[23,] -242.678808 -353.793021
[24,] 31.925883 -242.678808
[25,] -827.918197 31.925883
[26,] -776.363260 -827.918197
[27,] -604.625813 -776.363260
[28,] -221.849254 -604.625813
[29,] 107.942715 -221.849254
[30,] 211.103063 107.942715
[31,] 276.046542 211.103063
[32,] 354.279546 276.046542
[33,] 393.570691 354.279546
[34,] 450.984427 393.570691
[35,] 458.132638 450.984427
[36,] 717.054472 458.132638
[37,] 517.489065 717.054472
[38,] 339.134314 517.489065
[39,] 712.246663 339.134314
[40,] 814.005494 712.246663
[41,] 948.083462 814.005494
[42,] 902.320396 948.083462
[43,] 658.364677 902.320396
[44,] 891.038368 658.364677
[45,] 745.941065 891.038368
[46,] 578.762553 745.941065
[47,] 563.502211 578.762553
[48,] 426.747149 563.502211
[49,] 310.843359 426.747149
[50,] 331.259876 310.843359
[51,] 227.141189 331.259876
[52,] -7.222423 227.141189
[53,] 526.134572 -7.222423
[54,] 491.877270 526.134572
[55,] 376.866176 491.877270
[56,] 304.328662 376.866176
[57,] 127.633840 304.328662
[58,] -126.225291 127.633840
[59,] -344.543352 -126.225291
[60,] -454.120932 -344.543352
[61,] -535.620684 -454.120932
[62,] -604.468471 -535.620684
[63,] -689.490467 -604.468471
[64,] -727.072049 -689.490467
[65,] -832.253560 -727.072049
[66,] -756.626624 -832.253560
[67,] -757.040172 -756.626624
[68,] -770.451409 -757.040172
[69,] -851.554474 -770.451409
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 38.985168 51.474372
2 -164.230340 38.985168
3 -191.821336 -164.230340
4 171.136347 -191.821336
5 -77.772400 171.136347
6 -304.867312 -77.772400
7 -468.343303 -304.867312
8 -687.544262 -468.343303
9 -682.501773 -687.544262
10 -727.023147 -682.501773
11 -832.880037 -727.023147
12 480.702099 -832.880037
13 344.976209 480.702099
14 135.523871 344.976209
15 -73.534298 135.523871
16 -31.453916 -73.534298
17 38.773604 -31.453916
18 -139.661461 38.773604
19 -144.322497 -139.661461
20 118.374573 -144.322497
21 -164.832238 118.374573
22 -353.793021 -164.832238
23 -242.678808 -353.793021
24 31.925883 -242.678808
25 -827.918197 31.925883
26 -776.363260 -827.918197
27 -604.625813 -776.363260
28 -221.849254 -604.625813
29 107.942715 -221.849254
30 211.103063 107.942715
31 276.046542 211.103063
32 354.279546 276.046542
33 393.570691 354.279546
34 450.984427 393.570691
35 458.132638 450.984427
36 717.054472 458.132638
37 517.489065 717.054472
38 339.134314 517.489065
39 712.246663 339.134314
40 814.005494 712.246663
41 948.083462 814.005494
42 902.320396 948.083462
43 658.364677 902.320396
44 891.038368 658.364677
45 745.941065 891.038368
46 578.762553 745.941065
47 563.502211 578.762553
48 426.747149 563.502211
49 310.843359 426.747149
50 331.259876 310.843359
51 227.141189 331.259876
52 -7.222423 227.141189
53 526.134572 -7.222423
54 491.877270 526.134572
55 376.866176 491.877270
56 304.328662 376.866176
57 127.633840 304.328662
58 -126.225291 127.633840
59 -344.543352 -126.225291
60 -454.120932 -344.543352
61 -535.620684 -454.120932
62 -604.468471 -535.620684
63 -689.490467 -604.468471
64 -727.072049 -689.490467
65 -832.253560 -727.072049
66 -756.626624 -832.253560
67 -757.040172 -756.626624
68 -770.451409 -757.040172
69 -851.554474 -770.451409
> 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/freestat/rcomp/tmp/7e1ik1291108641.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/html/freestat/rcomp/tmp/8osin1291108641.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/html/freestat/rcomp/tmp/9osin1291108641.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/html/freestat/rcomp/tmp/10osin1291108641.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/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11k2xe1291108641.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/freestat/rcomp/tmp/1263wk1291108641.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/freestat/rcomp/tmp/13kdub1291108641.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/freestat/rcomp/tmp/14ndsy1291108641.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/freestat/rcomp/tmp/15g4911291108641.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/freestat/rcomp/tmp/16uepa1291108641.tab")
+ }
>
> try(system("convert tmp/1aj2w1291108641.ps tmp/1aj2w1291108641.png",intern=TRUE))
character(0)
> try(system("convert tmp/2aj2w1291108641.ps tmp/2aj2w1291108641.png",intern=TRUE))
character(0)
> try(system("convert tmp/3aj2w1291108641.ps tmp/3aj2w1291108641.png",intern=TRUE))
character(0)
> try(system("convert tmp/4la1z1291108641.ps tmp/4la1z1291108641.png",intern=TRUE))
character(0)
> try(system("convert tmp/5la1z1291108641.ps tmp/5la1z1291108641.png",intern=TRUE))
character(0)
> try(system("convert tmp/6la1z1291108641.ps tmp/6la1z1291108641.png",intern=TRUE))
character(0)
> try(system("convert tmp/7e1ik1291108641.ps tmp/7e1ik1291108641.png",intern=TRUE))
character(0)
> try(system("convert tmp/8osin1291108641.ps tmp/8osin1291108641.png",intern=TRUE))
character(0)
> try(system("convert tmp/9osin1291108641.ps tmp/9osin1291108641.png",intern=TRUE))
character(0)
> try(system("convert tmp/10osin1291108641.ps tmp/10osin1291108641.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.953 2.448 4.294