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(2756.76
+ ,10001.60
+ ,2849.27
+ ,10411.75
+ ,2921.44
+ ,10673.38
+ ,2981.85
+ ,10539.51
+ ,3080.58
+ ,10723.78
+ ,3106.22
+ ,10682.06
+ ,3119.31
+ ,10283.19
+ ,3061.26
+ ,10377.18
+ ,3097.31
+ ,10486.64
+ ,3161.69
+ ,10545.38
+ ,3257.16
+ ,10554.27
+ ,3277.01
+ ,10532.54
+ ,3295.32
+ ,10324.31
+ ,3363.99
+ ,10695.25
+ ,3494.17
+ ,10827.81
+ ,3667.03
+ ,10872.48
+ ,3813.06
+ ,10971.19
+ ,3917.96
+ ,11145.65
+ ,3895.51
+ ,11234.68
+ ,3801.06
+ ,11333.88
+ ,3570.12
+ ,10997.97
+ ,3701.61
+ ,11036.89
+ ,3862.27
+ ,11257.35
+ ,3970.10
+ ,11533.59
+ ,4138.52
+ ,11963.12
+ ,4199.75
+ ,12185.15
+ ,4290.89
+ ,12377.62
+ ,4443.91
+ ,12512.89
+ ,4502.64
+ ,12631.48
+ ,4356.98
+ ,12268.53
+ ,4591.27
+ ,12754.80
+ ,4696.96
+ ,13407.75
+ ,4621.40
+ ,13480.21
+ ,4562.84
+ ,13673.28
+ ,4202.52
+ ,13239.71
+ ,4296.49
+ ,13557.69
+ ,4435.23
+ ,13901.28
+ ,4105.18
+ ,13200.58
+ ,4116.68
+ ,13406.97
+ ,3844.49
+ ,12538.12
+ ,3720.98
+ ,12419.57
+ ,3674.40
+ ,12193.88
+ ,3857.62
+ ,12656.63
+ ,3801.06
+ ,12812.48
+ ,3504.37
+ ,12056.67
+ ,3032.60
+ ,11322.38
+ ,3047.03
+ ,11530.75
+ ,2962.34
+ ,11114.08
+ ,2197.82
+ ,9181.73
+ ,2014.45
+ ,8614.55
+ ,1862.83
+ ,8595.56
+ ,1905.41
+ ,8396.20
+ ,1810.99
+ ,7690.50
+ ,1670.07
+ ,7235.47
+ ,1864.44
+ ,7992.12
+ ,2052.02
+ ,8398.37
+ ,2029.60
+ ,8593.01
+ ,2070.83
+ ,8679.75
+ ,2293.41
+ ,9374.63
+ ,2443.27
+ ,9634.97)
+ ,dim=c(2
+ ,60)
+ ,dimnames=list(c('Bel20'
+ ,'Dow
')
+ ,1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Bel20','Dow
'),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 = 'Include Monthly Dummies'
> par1 = '2'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> 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
Dow\r Bel20 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 10001.60 2756.76 1 0 0 0 0 0 0 0 0 0 0
2 10411.75 2849.27 0 1 0 0 0 0 0 0 0 0 0
3 10673.38 2921.44 0 0 1 0 0 0 0 0 0 0 0
4 10539.51 2981.85 0 0 0 1 0 0 0 0 0 0 0
5 10723.78 3080.58 0 0 0 0 1 0 0 0 0 0 0
6 10682.06 3106.22 0 0 0 0 0 1 0 0 0 0 0
7 10283.19 3119.31 0 0 0 0 0 0 1 0 0 0 0
8 10377.18 3061.26 0 0 0 0 0 0 0 1 0 0 0
9 10486.64 3097.31 0 0 0 0 0 0 0 0 1 0 0
10 10545.38 3161.69 0 0 0 0 0 0 0 0 0 1 0
11 10554.27 3257.16 0 0 0 0 0 0 0 0 0 0 1
12 10532.54 3277.01 0 0 0 0 0 0 0 0 0 0 0
13 10324.31 3295.32 1 0 0 0 0 0 0 0 0 0 0
14 10695.25 3363.99 0 1 0 0 0 0 0 0 0 0 0
15 10827.81 3494.17 0 0 1 0 0 0 0 0 0 0 0
16 10872.48 3667.03 0 0 0 1 0 0 0 0 0 0 0
17 10971.19 3813.06 0 0 0 0 1 0 0 0 0 0 0
18 11145.65 3917.96 0 0 0 0 0 1 0 0 0 0 0
19 11234.68 3895.51 0 0 0 0 0 0 1 0 0 0 0
20 11333.88 3801.06 0 0 0 0 0 0 0 1 0 0 0
21 10997.97 3570.12 0 0 0 0 0 0 0 0 1 0 0
22 11036.89 3701.61 0 0 0 0 0 0 0 0 0 1 0
23 11257.35 3862.27 0 0 0 0 0 0 0 0 0 0 1
24 11533.59 3970.10 0 0 0 0 0 0 0 0 0 0 0
25 11963.12 4138.52 1 0 0 0 0 0 0 0 0 0 0
26 12185.15 4199.75 0 1 0 0 0 0 0 0 0 0 0
27 12377.62 4290.89 0 0 1 0 0 0 0 0 0 0 0
28 12512.89 4443.91 0 0 0 1 0 0 0 0 0 0 0
29 12631.48 4502.64 0 0 0 0 1 0 0 0 0 0 0
30 12268.53 4356.98 0 0 0 0 0 1 0 0 0 0 0
31 12754.80 4591.27 0 0 0 0 0 0 1 0 0 0 0
32 13407.75 4696.96 0 0 0 0 0 0 0 1 0 0 0
33 13480.21 4621.40 0 0 0 0 0 0 0 0 1 0 0
34 13673.28 4562.84 0 0 0 0 0 0 0 0 0 1 0
35 13239.71 4202.52 0 0 0 0 0 0 0 0 0 0 1
36 13557.69 4296.49 0 0 0 0 0 0 0 0 0 0 0
37 13901.28 4435.23 1 0 0 0 0 0 0 0 0 0 0
38 13200.58 4105.18 0 1 0 0 0 0 0 0 0 0 0
39 13406.97 4116.68 0 0 1 0 0 0 0 0 0 0 0
40 12538.12 3844.49 0 0 0 1 0 0 0 0 0 0 0
41 12419.57 3720.98 0 0 0 0 1 0 0 0 0 0 0
42 12193.88 3674.40 0 0 0 0 0 1 0 0 0 0 0
43 12656.63 3857.62 0 0 0 0 0 0 1 0 0 0 0
44 12812.48 3801.06 0 0 0 0 0 0 0 1 0 0 0
45 12056.67 3504.37 0 0 0 0 0 0 0 0 1 0 0
46 11322.38 3032.60 0 0 0 0 0 0 0 0 0 1 0
47 11530.75 3047.03 0 0 0 0 0 0 0 0 0 0 1
48 11114.08 2962.34 0 0 0 0 0 0 0 0 0 0 0
49 9181.73 2197.82 1 0 0 0 0 0 0 0 0 0 0
50 8614.55 2014.45 0 1 0 0 0 0 0 0 0 0 0
51 8595.56 1862.83 0 0 1 0 0 0 0 0 0 0 0
52 8396.20 1905.41 0 0 0 1 0 0 0 0 0 0 0
53 7690.50 1810.99 0 0 0 0 1 0 0 0 0 0 0
54 7235.47 1670.07 0 0 0 0 0 1 0 0 0 0 0
55 7992.12 1864.44 0 0 0 0 0 0 1 0 0 0 0
56 8398.37 2052.02 0 0 0 0 0 0 0 1 0 0 0
57 8593.01 2029.60 0 0 0 0 0 0 0 0 1 0 0
58 8679.75 2070.83 0 0 0 0 0 0 0 0 0 1 0
59 9374.63 2293.41 0 0 0 0 0 0 0 0 0 0 1
60 9634.97 2443.27 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Bel20 M1 M2 M3 M4
5080.608 1.827 -154.281 -100.886 -2.121 -263.807
M5 M6 M7 M8 M9 M10
-379.610 -487.750 -428.771 -177.897 -105.478 -69.684
M11
21.584
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-902.0 -476.2 -72.8 398.8 964.4
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5080.60794 408.72732 12.430 <2e-16 ***
Bel20 1.82721 0.09147 19.977 <2e-16 ***
M1 -154.28101 376.62769 -0.410 0.684
M2 -100.88551 376.69778 -0.268 0.790
M3 -2.12147 376.65147 -0.006 0.996
M4 -263.80704 376.62573 -0.700 0.487
M5 -379.61032 376.62089 -1.008 0.319
M6 -487.75031 376.64290 -1.295 0.202
M7 -428.77088 376.68448 -1.138 0.261
M8 -177.89681 376.71598 -0.472 0.639
M9 -105.47838 376.62779 -0.280 0.781
M10 -69.68360 376.69892 -0.185 0.854
M11 21.58429 376.65724 0.057 0.955
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 595.5 on 47 degrees of freedom
Multiple R-squared: 0.8956, Adjusted R-squared: 0.869
F-statistic: 33.6 on 12 and 47 DF, p-value: < 2.2e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 9.152971e-04 1.830594e-03 0.9990847029
[2,] 1.212794e-04 2.425588e-04 0.9998787206
[3,] 1.977907e-05 3.955814e-05 0.9999802209
[4,] 9.825212e-04 1.965042e-03 0.9990174788
[5,] 1.273467e-03 2.546935e-03 0.9987265326
[6,] 4.460221e-04 8.920443e-04 0.9995539779
[7,] 1.631172e-04 3.262344e-04 0.9998368828
[8,] 1.223053e-04 2.446106e-04 0.9998776947
[9,] 2.421079e-04 4.842159e-04 0.9997578921
[10,] 3.249028e-03 6.498055e-03 0.9967509724
[11,] 4.332594e-03 8.665189e-03 0.9956674055
[12,] 6.480965e-03 1.296193e-02 0.9935190355
[13,] 1.144350e-02 2.288700e-02 0.9885564979
[14,] 1.780061e-02 3.560122e-02 0.9821993907
[15,] 2.108661e-02 4.217322e-02 0.9789133884
[16,] 7.657437e-02 1.531487e-01 0.9234256345
[17,] 2.527450e-01 5.054899e-01 0.7472550485
[18,] 5.404964e-01 9.190072e-01 0.4595035916
[19,] 8.529608e-01 2.940784e-01 0.1470392194
[20,] 9.699210e-01 6.015803e-02 0.0300790159
[21,] 9.902462e-01 1.950758e-02 0.0097537904
[22,] 9.954049e-01 9.190186e-03 0.0045950928
[23,] 9.945224e-01 1.095530e-02 0.0054776491
[24,] 9.975669e-01 4.866239e-03 0.0024331194
[25,] 9.994946e-01 1.010792e-03 0.0005053962
[26,] 9.984073e-01 3.185456e-03 0.0015927279
[27,] 9.947050e-01 1.059005e-02 0.0052950233
[28,] 9.921946e-01 1.561090e-02 0.0078054483
[29,] 9.705385e-01 5.892297e-02 0.0294614847
> postscript(file="/var/www/html/rcomp/tmp/152o71259618846.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/2p9bj1259618846.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/3fyah1259618846.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/4bfvo1259618846.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/5rgwi1259618846.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 = 60
Frequency = 1
1 2 3 4 5 6 7
38.08317 225.80212 256.79806 274.23165 393.90411 413.47434 -68.29332
8 9 10 11 12 13 14
-119.10762 -147.93711 -242.62792 -499.44991 -535.86581 -623.27109 -431.20136
15 16 17 18 19 20 21
-635.27209 -644.76870 -697.08345 -606.15819 -535.08667 -514.18039 -500.53207
22 23 24 25 26 27 28
-737.66719 -902.03525 -801.23942 -525.16776 -468.41356 -541.23987 -423.88455
29 30 31 32 33 34 35
-296.80353 -285.46159 -286.26893 -77.31122 60.79462 325.07148 458.61526
36 37 38 39 40 41 42
626.47627 870.83963 719.81605 806.42905 696.61394 919.54640 887.10801
43 44 45 46 47 48 49
956.09647 964.41961 678.30724 770.24711 860.98252 620.64355 239.51605
50 51 52 53 54 55 56
-46.00326 113.28485 97.80766 -319.56353 -408.96257 -66.44755 -253.82038
57 58 59 60
-90.63268 -115.02348 81.88738 89.98541
> postscript(file="/var/www/html/rcomp/tmp/66r711259618846.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 38.08317 NA
1 225.80212 38.08317
2 256.79806 225.80212
3 274.23165 256.79806
4 393.90411 274.23165
5 413.47434 393.90411
6 -68.29332 413.47434
7 -119.10762 -68.29332
8 -147.93711 -119.10762
9 -242.62792 -147.93711
10 -499.44991 -242.62792
11 -535.86581 -499.44991
12 -623.27109 -535.86581
13 -431.20136 -623.27109
14 -635.27209 -431.20136
15 -644.76870 -635.27209
16 -697.08345 -644.76870
17 -606.15819 -697.08345
18 -535.08667 -606.15819
19 -514.18039 -535.08667
20 -500.53207 -514.18039
21 -737.66719 -500.53207
22 -902.03525 -737.66719
23 -801.23942 -902.03525
24 -525.16776 -801.23942
25 -468.41356 -525.16776
26 -541.23987 -468.41356
27 -423.88455 -541.23987
28 -296.80353 -423.88455
29 -285.46159 -296.80353
30 -286.26893 -285.46159
31 -77.31122 -286.26893
32 60.79462 -77.31122
33 325.07148 60.79462
34 458.61526 325.07148
35 626.47627 458.61526
36 870.83963 626.47627
37 719.81605 870.83963
38 806.42905 719.81605
39 696.61394 806.42905
40 919.54640 696.61394
41 887.10801 919.54640
42 956.09647 887.10801
43 964.41961 956.09647
44 678.30724 964.41961
45 770.24711 678.30724
46 860.98252 770.24711
47 620.64355 860.98252
48 239.51605 620.64355
49 -46.00326 239.51605
50 113.28485 -46.00326
51 97.80766 113.28485
52 -319.56353 97.80766
53 -408.96257 -319.56353
54 -66.44755 -408.96257
55 -253.82038 -66.44755
56 -90.63268 -253.82038
57 -115.02348 -90.63268
58 81.88738 -115.02348
59 89.98541 81.88738
60 NA 89.98541
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 225.80212 38.08317
[2,] 256.79806 225.80212
[3,] 274.23165 256.79806
[4,] 393.90411 274.23165
[5,] 413.47434 393.90411
[6,] -68.29332 413.47434
[7,] -119.10762 -68.29332
[8,] -147.93711 -119.10762
[9,] -242.62792 -147.93711
[10,] -499.44991 -242.62792
[11,] -535.86581 -499.44991
[12,] -623.27109 -535.86581
[13,] -431.20136 -623.27109
[14,] -635.27209 -431.20136
[15,] -644.76870 -635.27209
[16,] -697.08345 -644.76870
[17,] -606.15819 -697.08345
[18,] -535.08667 -606.15819
[19,] -514.18039 -535.08667
[20,] -500.53207 -514.18039
[21,] -737.66719 -500.53207
[22,] -902.03525 -737.66719
[23,] -801.23942 -902.03525
[24,] -525.16776 -801.23942
[25,] -468.41356 -525.16776
[26,] -541.23987 -468.41356
[27,] -423.88455 -541.23987
[28,] -296.80353 -423.88455
[29,] -285.46159 -296.80353
[30,] -286.26893 -285.46159
[31,] -77.31122 -286.26893
[32,] 60.79462 -77.31122
[33,] 325.07148 60.79462
[34,] 458.61526 325.07148
[35,] 626.47627 458.61526
[36,] 870.83963 626.47627
[37,] 719.81605 870.83963
[38,] 806.42905 719.81605
[39,] 696.61394 806.42905
[40,] 919.54640 696.61394
[41,] 887.10801 919.54640
[42,] 956.09647 887.10801
[43,] 964.41961 956.09647
[44,] 678.30724 964.41961
[45,] 770.24711 678.30724
[46,] 860.98252 770.24711
[47,] 620.64355 860.98252
[48,] 239.51605 620.64355
[49,] -46.00326 239.51605
[50,] 113.28485 -46.00326
[51,] 97.80766 113.28485
[52,] -319.56353 97.80766
[53,] -408.96257 -319.56353
[54,] -66.44755 -408.96257
[55,] -253.82038 -66.44755
[56,] -90.63268 -253.82038
[57,] -115.02348 -90.63268
[58,] 81.88738 -115.02348
[59,] 89.98541 81.88738
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 225.80212 38.08317
2 256.79806 225.80212
3 274.23165 256.79806
4 393.90411 274.23165
5 413.47434 393.90411
6 -68.29332 413.47434
7 -119.10762 -68.29332
8 -147.93711 -119.10762
9 -242.62792 -147.93711
10 -499.44991 -242.62792
11 -535.86581 -499.44991
12 -623.27109 -535.86581
13 -431.20136 -623.27109
14 -635.27209 -431.20136
15 -644.76870 -635.27209
16 -697.08345 -644.76870
17 -606.15819 -697.08345
18 -535.08667 -606.15819
19 -514.18039 -535.08667
20 -500.53207 -514.18039
21 -737.66719 -500.53207
22 -902.03525 -737.66719
23 -801.23942 -902.03525
24 -525.16776 -801.23942
25 -468.41356 -525.16776
26 -541.23987 -468.41356
27 -423.88455 -541.23987
28 -296.80353 -423.88455
29 -285.46159 -296.80353
30 -286.26893 -285.46159
31 -77.31122 -286.26893
32 60.79462 -77.31122
33 325.07148 60.79462
34 458.61526 325.07148
35 626.47627 458.61526
36 870.83963 626.47627
37 719.81605 870.83963
38 806.42905 719.81605
39 696.61394 806.42905
40 919.54640 696.61394
41 887.10801 919.54640
42 956.09647 887.10801
43 964.41961 956.09647
44 678.30724 964.41961
45 770.24711 678.30724
46 860.98252 770.24711
47 620.64355 860.98252
48 239.51605 620.64355
49 -46.00326 239.51605
50 113.28485 -46.00326
51 97.80766 113.28485
52 -319.56353 97.80766
53 -408.96257 -319.56353
54 -66.44755 -408.96257
55 -253.82038 -66.44755
56 -90.63268 -253.82038
57 -115.02348 -90.63268
58 81.88738 -115.02348
59 89.98541 81.88738
> 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/7c6251259618846.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/8k9321259618846.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/9v0b31259618846.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/10kppn1259618846.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/11ze3y1259618846.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/12rqei1259618846.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/135szf1259618846.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/14h70k1259618846.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/159q0h1259618846.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/16axgv1259618846.tab")
+ }
>
> system("convert tmp/152o71259618846.ps tmp/152o71259618846.png")
> system("convert tmp/2p9bj1259618846.ps tmp/2p9bj1259618846.png")
> system("convert tmp/3fyah1259618846.ps tmp/3fyah1259618846.png")
> system("convert tmp/4bfvo1259618846.ps tmp/4bfvo1259618846.png")
> system("convert tmp/5rgwi1259618846.ps tmp/5rgwi1259618846.png")
> system("convert tmp/66r711259618846.ps tmp/66r711259618846.png")
> system("convert tmp/7c6251259618846.ps tmp/7c6251259618846.png")
> system("convert tmp/8k9321259618846.ps tmp/8k9321259618846.png")
> system("convert tmp/9v0b31259618846.ps tmp/9v0b31259618846.png")
> system("convert tmp/10kppn1259618846.ps tmp/10kppn1259618846.png")
>
>
> proc.time()
user system elapsed
2.378 1.531 3.176