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.
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(15044.5,1,14944.2,1,16754.8,1,14254,1,15454.9,1,15644.8,1,14568.3,1,12520.2,1,14803,1,15873.2,1,14755.3,1,12875.1,1,14291.1,1,14205.3,1,15859.4,1,15258.9,1,15498.6,1,15106.5,1,15023.6,1,12083,1,15761.3,1,16943,1,15070.3,1,13659.6,1,14768.9,0,14725.1,0,15998.1,0,15370.6,0,14956.9,0,15469.7,0,15101.8,0,11703.7,0,16283.6,0,16726.5,0,14968.9,0,14861,0,14583.3,0,15305.8,0,17903.9,0,16379.4,0,15420.3,0,17870.5,0,15912.8,0,13866.5,0,17823.2,0,17872,0,17420.4,0,16704.4,0,15991.2,0,16583.6,0,19123.5,0,17838.7,0,17209.4,0,18586.5,0,16258.1,0,15141.6,0,19202.1,0,17746.5,0,19090.1,0,18040.3,0,17515.5,0,17751.8,0,21072.4,0,17170,0,19439.5,0,19795.4,0,17574.9,0,16165.4,0,19464.6,0,19932.1,0,19961.2,0,17343.4,0,18924.2,0,18574.1,0,21350.6,0,18594.6,0,19823.1,0,20844.4,0,19640.2,0,17735.4,0,19813.6,0,22160,0,20664.3,0,17877.4,0,21211.2,0,21423.1,0,21688.7,0,23243.2,0,21490.2,0,22925.8,0,23184.8,0,18562.2,0),dim=c(2,92),dimnames=list(c('Y','X'),1:92))
> y <- array(NA,dim=c(2,92),dimnames=list(c('Y','X'),1:92))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par20 = ''
> par19 = ''
> par18 = ''
> par17 = ''
> par16 = ''
> par15 = ''
> par14 = ''
> par13 = ''
> par12 = ''
> par11 = ''
> par10 = ''
> par9 = ''
> par8 = ''
> par7 = ''
> par6 = ''
> par5 = ''
> par4 = ''
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> ylab = ''
> xlab = ''
> main = ''
> #'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
Y X
1 15044.5 1
2 14944.2 1
3 16754.8 1
4 14254.0 1
5 15454.9 1
6 15644.8 1
7 14568.3 1
8 12520.2 1
9 14803.0 1
10 15873.2 1
11 14755.3 1
12 12875.1 1
13 14291.1 1
14 14205.3 1
15 15859.4 1
16 15258.9 1
17 15498.6 1
18 15106.5 1
19 15023.6 1
20 12083.0 1
21 15761.3 1
22 16943.0 1
23 15070.3 1
24 13659.6 1
25 14768.9 0
26 14725.1 0
27 15998.1 0
28 15370.6 0
29 14956.9 0
30 15469.7 0
31 15101.8 0
32 11703.7 0
33 16283.6 0
34 16726.5 0
35 14968.9 0
36 14861.0 0
37 14583.3 0
38 15305.8 0
39 17903.9 0
40 16379.4 0
41 15420.3 0
42 17870.5 0
43 15912.8 0
44 13866.5 0
45 17823.2 0
46 17872.0 0
47 17420.4 0
48 16704.4 0
49 15991.2 0
50 16583.6 0
51 19123.5 0
52 17838.7 0
53 17209.4 0
54 18586.5 0
55 16258.1 0
56 15141.6 0
57 19202.1 0
58 17746.5 0
59 19090.1 0
60 18040.3 0
61 17515.5 0
62 17751.8 0
63 21072.4 0
64 17170.0 0
65 19439.5 0
66 19795.4 0
67 17574.9 0
68 16165.4 0
69 19464.6 0
70 19932.1 0
71 19961.2 0
72 17343.4 0
73 18924.2 0
74 18574.1 0
75 21350.6 0
76 18594.6 0
77 19823.1 0
78 20844.4 0
79 19640.2 0
80 17735.4 0
81 19813.6 0
82 22160.0 0
83 20664.3 0
84 17877.4 0
85 21211.2 0
86 21423.1 0
87 21688.7 0
88 23243.2 0
89 21490.2 0
90 22925.8 0
91 23184.8 0
92 18562.2 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X
17967 -3123
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-6262.86 -1611.11 -88.87 1176.59 5276.64
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17966.6 267.0 67.284 < 2e-16 ***
X -3122.7 522.8 -5.973 4.57e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2202 on 90 degrees of freedom
Multiple R-squared: 0.2839, Adjusted R-squared: 0.2759
F-statistic: 35.68 on 1 and 90 DF, p-value: 4.575e-08
> 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,] 1.211326e-01 2.422653e-01 0.8788674
[2,] 4.790053e-02 9.580106e-02 0.9520995
[3,] 2.241747e-02 4.483495e-02 0.9775825
[4,] 8.485218e-02 1.697044e-01 0.9151478
[5,] 4.229545e-02 8.459090e-02 0.9577045
[6,] 2.565868e-02 5.131736e-02 0.9743413
[7,] 1.197317e-02 2.394633e-02 0.9880268
[8,] 1.752620e-02 3.505239e-02 0.9824738
[9,] 9.036742e-03 1.807348e-02 0.9909633
[10,] 4.591251e-03 9.182502e-03 0.9954087
[11,] 3.036614e-03 6.073228e-03 0.9969634
[12,] 1.455416e-03 2.910832e-03 0.9985446
[13,] 7.360504e-04 1.472101e-03 0.9992639
[14,] 3.188925e-04 6.377850e-04 0.9996811
[15,] 1.319598e-04 2.639195e-04 0.9998680
[16,] 7.108197e-04 1.421639e-03 0.9992892
[17,] 4.315821e-04 8.631642e-04 0.9995684
[18,] 6.741315e-04 1.348263e-03 0.9993259
[19,] 3.303983e-04 6.607966e-04 0.9996696
[20,] 2.237956e-04 4.475913e-04 0.9997762
[21,] 1.289575e-04 2.579150e-04 0.9998710
[22,] 7.553795e-05 1.510759e-04 0.9999245
[23,] 5.072731e-05 1.014546e-04 0.9999493
[24,] 2.780027e-05 5.560054e-05 0.9999722
[25,] 1.651211e-05 3.302422e-05 0.9999835
[26,] 9.185584e-06 1.837117e-05 0.9999908
[27,] 5.408535e-06 1.081707e-05 0.9999946
[28,] 1.874878e-04 3.749756e-04 0.9998125
[29,] 1.760015e-04 3.520029e-04 0.9998240
[30,] 1.864543e-04 3.729086e-04 0.9998135
[31,] 1.522899e-04 3.045798e-04 0.9998477
[32,] 1.356863e-04 2.713727e-04 0.9998643
[33,] 1.468532e-04 2.937064e-04 0.9998531
[34,] 1.306245e-04 2.612490e-04 0.9998694
[35,] 3.350756e-04 6.701513e-04 0.9996649
[36,] 2.987718e-04 5.975436e-04 0.9997012
[37,] 2.874791e-04 5.749582e-04 0.9997125
[38,] 4.969368e-04 9.938735e-04 0.9995031
[39,] 4.557517e-04 9.115035e-04 0.9995442
[40,] 1.587604e-03 3.175209e-03 0.9984124
[41,] 2.307268e-03 4.614536e-03 0.9976927
[42,] 3.049202e-03 6.098404e-03 0.9969508
[43,] 3.222608e-03 6.445216e-03 0.9967774
[44,] 3.148025e-03 6.296049e-03 0.9968520
[45,] 3.639977e-03 7.279954e-03 0.9963600
[46,] 3.873997e-03 7.747994e-03 0.9961260
[47,] 7.766044e-03 1.553209e-02 0.9922340
[48,] 8.147659e-03 1.629532e-02 0.9918523
[49,] 8.122175e-03 1.624435e-02 0.9918778
[50,] 9.746246e-03 1.949249e-02 0.9902538
[51,] 1.226513e-02 2.453026e-02 0.9877349
[52,] 3.130558e-02 6.261116e-02 0.9686944
[53,] 4.180366e-02 8.360733e-02 0.9581963
[54,] 4.350549e-02 8.701098e-02 0.9564945
[55,] 5.037065e-02 1.007413e-01 0.9496294
[56,] 5.078986e-02 1.015797e-01 0.9492101
[57,] 5.459188e-02 1.091838e-01 0.9454081
[58,] 5.797476e-02 1.159495e-01 0.9420252
[59,] 1.123940e-01 2.247881e-01 0.8876060
[60,] 1.324926e-01 2.649851e-01 0.8675074
[61,] 1.340196e-01 2.680391e-01 0.8659804
[62,] 1.378917e-01 2.757835e-01 0.8621083
[63,] 1.515573e-01 3.031145e-01 0.8484427
[64,] 2.875175e-01 5.750350e-01 0.7124825
[65,] 2.796640e-01 5.593280e-01 0.7203360
[66,] 2.732101e-01 5.464201e-01 0.7267899
[67,] 2.614091e-01 5.228182e-01 0.7385909
[68,] 3.431858e-01 6.863716e-01 0.6568142
[69,] 3.380406e-01 6.760811e-01 0.6619594
[70,] 3.555393e-01 7.110787e-01 0.6444607
[71,] 3.734020e-01 7.468039e-01 0.6265980
[72,] 3.920205e-01 7.840409e-01 0.6079795
[73,] 3.599885e-01 7.199769e-01 0.6400115
[74,] 3.283726e-01 6.567452e-01 0.6716274
[75,] 2.957559e-01 5.915117e-01 0.7042441
[76,] 4.459091e-01 8.918182e-01 0.5540909
[77,] 4.201826e-01 8.403651e-01 0.5798174
[78,] 4.030957e-01 8.061914e-01 0.5969043
[79,] 3.332422e-01 6.664844e-01 0.6667578
[80,] 6.026531e-01 7.946938e-01 0.3973469
[81,] 5.069987e-01 9.860026e-01 0.4930013
[82,] 3.938306e-01 7.876613e-01 0.6061694
[83,] 2.712247e-01 5.424495e-01 0.7287753
> postscript(file="/var/www/html/rcomp/tmp/19ulz1229016224.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/2fea61229016224.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/39m7w1229016224.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/4xime1229016224.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/5v1kj1229016224.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 = 92
Frequency = 1
1 2 3 4 5 6
200.62917 100.32917 1910.92917 -589.87083 611.02917 800.92917
7 8 9 10 11 12
-275.57083 -2323.67083 -40.87083 1029.32917 -88.57083 -1968.77083
13 14 15 16 17 18
-552.77083 -638.57083 1015.52917 415.02917 654.72917 262.62917
19 20 21 22 23 24
179.72917 -2760.87083 917.42917 2099.12917 226.42917 -1184.27083
25 26 27 28 29 30
-3197.66176 -3241.46176 -1968.46176 -2595.96176 -3009.66176 -2496.86176
31 32 33 34 35 36
-2864.76176 -6262.86176 -1682.96176 -1240.06176 -2997.66176 -3105.56176
37 38 39 40 41 42
-3383.26176 -2660.76176 -62.66176 -1587.16176 -2546.26176 -96.06176
43 44 45 46 47 48
-2053.76176 -4100.06176 -143.36176 -94.56176 -546.16176 -1262.16176
49 50 51 52 53 54
-1975.36176 -1382.96176 1156.93824 -127.86176 -757.16176 619.93824
55 56 57 58 59 60
-1708.46176 -2824.96176 1235.53824 -220.06176 1123.53824 73.73824
61 62 63 64 65 66
-451.06176 -214.76176 3105.83824 -796.56176 1472.93824 1828.83824
67 68 69 70 71 72
-391.66176 -1801.16176 1498.03824 1965.53824 1994.63824 -623.16176
73 74 75 76 77 78
957.63824 607.53824 3384.03824 628.03824 1856.53824 2877.83824
79 80 81 82 83 84
1673.63824 -231.16176 1847.03824 4193.43824 2697.73824 -89.16176
85 86 87 88 89 90
3244.63824 3456.53824 3722.13824 5276.63824 3523.63824 4959.23824
91 92
5218.23824 595.63824
> postscript(file="/var/www/html/rcomp/tmp/6vdrs1229016224.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 = 92
Frequency = 1
lag(myerror, k = 1) myerror
0 200.62917 NA
1 100.32917 200.62917
2 1910.92917 100.32917
3 -589.87083 1910.92917
4 611.02917 -589.87083
5 800.92917 611.02917
6 -275.57083 800.92917
7 -2323.67083 -275.57083
8 -40.87083 -2323.67083
9 1029.32917 -40.87083
10 -88.57083 1029.32917
11 -1968.77083 -88.57083
12 -552.77083 -1968.77083
13 -638.57083 -552.77083
14 1015.52917 -638.57083
15 415.02917 1015.52917
16 654.72917 415.02917
17 262.62917 654.72917
18 179.72917 262.62917
19 -2760.87083 179.72917
20 917.42917 -2760.87083
21 2099.12917 917.42917
22 226.42917 2099.12917
23 -1184.27083 226.42917
24 -3197.66176 -1184.27083
25 -3241.46176 -3197.66176
26 -1968.46176 -3241.46176
27 -2595.96176 -1968.46176
28 -3009.66176 -2595.96176
29 -2496.86176 -3009.66176
30 -2864.76176 -2496.86176
31 -6262.86176 -2864.76176
32 -1682.96176 -6262.86176
33 -1240.06176 -1682.96176
34 -2997.66176 -1240.06176
35 -3105.56176 -2997.66176
36 -3383.26176 -3105.56176
37 -2660.76176 -3383.26176
38 -62.66176 -2660.76176
39 -1587.16176 -62.66176
40 -2546.26176 -1587.16176
41 -96.06176 -2546.26176
42 -2053.76176 -96.06176
43 -4100.06176 -2053.76176
44 -143.36176 -4100.06176
45 -94.56176 -143.36176
46 -546.16176 -94.56176
47 -1262.16176 -546.16176
48 -1975.36176 -1262.16176
49 -1382.96176 -1975.36176
50 1156.93824 -1382.96176
51 -127.86176 1156.93824
52 -757.16176 -127.86176
53 619.93824 -757.16176
54 -1708.46176 619.93824
55 -2824.96176 -1708.46176
56 1235.53824 -2824.96176
57 -220.06176 1235.53824
58 1123.53824 -220.06176
59 73.73824 1123.53824
60 -451.06176 73.73824
61 -214.76176 -451.06176
62 3105.83824 -214.76176
63 -796.56176 3105.83824
64 1472.93824 -796.56176
65 1828.83824 1472.93824
66 -391.66176 1828.83824
67 -1801.16176 -391.66176
68 1498.03824 -1801.16176
69 1965.53824 1498.03824
70 1994.63824 1965.53824
71 -623.16176 1994.63824
72 957.63824 -623.16176
73 607.53824 957.63824
74 3384.03824 607.53824
75 628.03824 3384.03824
76 1856.53824 628.03824
77 2877.83824 1856.53824
78 1673.63824 2877.83824
79 -231.16176 1673.63824
80 1847.03824 -231.16176
81 4193.43824 1847.03824
82 2697.73824 4193.43824
83 -89.16176 2697.73824
84 3244.63824 -89.16176
85 3456.53824 3244.63824
86 3722.13824 3456.53824
87 5276.63824 3722.13824
88 3523.63824 5276.63824
89 4959.23824 3523.63824
90 5218.23824 4959.23824
91 595.63824 5218.23824
92 NA 595.63824
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 100.32917 200.62917
[2,] 1910.92917 100.32917
[3,] -589.87083 1910.92917
[4,] 611.02917 -589.87083
[5,] 800.92917 611.02917
[6,] -275.57083 800.92917
[7,] -2323.67083 -275.57083
[8,] -40.87083 -2323.67083
[9,] 1029.32917 -40.87083
[10,] -88.57083 1029.32917
[11,] -1968.77083 -88.57083
[12,] -552.77083 -1968.77083
[13,] -638.57083 -552.77083
[14,] 1015.52917 -638.57083
[15,] 415.02917 1015.52917
[16,] 654.72917 415.02917
[17,] 262.62917 654.72917
[18,] 179.72917 262.62917
[19,] -2760.87083 179.72917
[20,] 917.42917 -2760.87083
[21,] 2099.12917 917.42917
[22,] 226.42917 2099.12917
[23,] -1184.27083 226.42917
[24,] -3197.66176 -1184.27083
[25,] -3241.46176 -3197.66176
[26,] -1968.46176 -3241.46176
[27,] -2595.96176 -1968.46176
[28,] -3009.66176 -2595.96176
[29,] -2496.86176 -3009.66176
[30,] -2864.76176 -2496.86176
[31,] -6262.86176 -2864.76176
[32,] -1682.96176 -6262.86176
[33,] -1240.06176 -1682.96176
[34,] -2997.66176 -1240.06176
[35,] -3105.56176 -2997.66176
[36,] -3383.26176 -3105.56176
[37,] -2660.76176 -3383.26176
[38,] -62.66176 -2660.76176
[39,] -1587.16176 -62.66176
[40,] -2546.26176 -1587.16176
[41,] -96.06176 -2546.26176
[42,] -2053.76176 -96.06176
[43,] -4100.06176 -2053.76176
[44,] -143.36176 -4100.06176
[45,] -94.56176 -143.36176
[46,] -546.16176 -94.56176
[47,] -1262.16176 -546.16176
[48,] -1975.36176 -1262.16176
[49,] -1382.96176 -1975.36176
[50,] 1156.93824 -1382.96176
[51,] -127.86176 1156.93824
[52,] -757.16176 -127.86176
[53,] 619.93824 -757.16176
[54,] -1708.46176 619.93824
[55,] -2824.96176 -1708.46176
[56,] 1235.53824 -2824.96176
[57,] -220.06176 1235.53824
[58,] 1123.53824 -220.06176
[59,] 73.73824 1123.53824
[60,] -451.06176 73.73824
[61,] -214.76176 -451.06176
[62,] 3105.83824 -214.76176
[63,] -796.56176 3105.83824
[64,] 1472.93824 -796.56176
[65,] 1828.83824 1472.93824
[66,] -391.66176 1828.83824
[67,] -1801.16176 -391.66176
[68,] 1498.03824 -1801.16176
[69,] 1965.53824 1498.03824
[70,] 1994.63824 1965.53824
[71,] -623.16176 1994.63824
[72,] 957.63824 -623.16176
[73,] 607.53824 957.63824
[74,] 3384.03824 607.53824
[75,] 628.03824 3384.03824
[76,] 1856.53824 628.03824
[77,] 2877.83824 1856.53824
[78,] 1673.63824 2877.83824
[79,] -231.16176 1673.63824
[80,] 1847.03824 -231.16176
[81,] 4193.43824 1847.03824
[82,] 2697.73824 4193.43824
[83,] -89.16176 2697.73824
[84,] 3244.63824 -89.16176
[85,] 3456.53824 3244.63824
[86,] 3722.13824 3456.53824
[87,] 5276.63824 3722.13824
[88,] 3523.63824 5276.63824
[89,] 4959.23824 3523.63824
[90,] 5218.23824 4959.23824
[91,] 595.63824 5218.23824
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 100.32917 200.62917
2 1910.92917 100.32917
3 -589.87083 1910.92917
4 611.02917 -589.87083
5 800.92917 611.02917
6 -275.57083 800.92917
7 -2323.67083 -275.57083
8 -40.87083 -2323.67083
9 1029.32917 -40.87083
10 -88.57083 1029.32917
11 -1968.77083 -88.57083
12 -552.77083 -1968.77083
13 -638.57083 -552.77083
14 1015.52917 -638.57083
15 415.02917 1015.52917
16 654.72917 415.02917
17 262.62917 654.72917
18 179.72917 262.62917
19 -2760.87083 179.72917
20 917.42917 -2760.87083
21 2099.12917 917.42917
22 226.42917 2099.12917
23 -1184.27083 226.42917
24 -3197.66176 -1184.27083
25 -3241.46176 -3197.66176
26 -1968.46176 -3241.46176
27 -2595.96176 -1968.46176
28 -3009.66176 -2595.96176
29 -2496.86176 -3009.66176
30 -2864.76176 -2496.86176
31 -6262.86176 -2864.76176
32 -1682.96176 -6262.86176
33 -1240.06176 -1682.96176
34 -2997.66176 -1240.06176
35 -3105.56176 -2997.66176
36 -3383.26176 -3105.56176
37 -2660.76176 -3383.26176
38 -62.66176 -2660.76176
39 -1587.16176 -62.66176
40 -2546.26176 -1587.16176
41 -96.06176 -2546.26176
42 -2053.76176 -96.06176
43 -4100.06176 -2053.76176
44 -143.36176 -4100.06176
45 -94.56176 -143.36176
46 -546.16176 -94.56176
47 -1262.16176 -546.16176
48 -1975.36176 -1262.16176
49 -1382.96176 -1975.36176
50 1156.93824 -1382.96176
51 -127.86176 1156.93824
52 -757.16176 -127.86176
53 619.93824 -757.16176
54 -1708.46176 619.93824
55 -2824.96176 -1708.46176
56 1235.53824 -2824.96176
57 -220.06176 1235.53824
58 1123.53824 -220.06176
59 73.73824 1123.53824
60 -451.06176 73.73824
61 -214.76176 -451.06176
62 3105.83824 -214.76176
63 -796.56176 3105.83824
64 1472.93824 -796.56176
65 1828.83824 1472.93824
66 -391.66176 1828.83824
67 -1801.16176 -391.66176
68 1498.03824 -1801.16176
69 1965.53824 1498.03824
70 1994.63824 1965.53824
71 -623.16176 1994.63824
72 957.63824 -623.16176
73 607.53824 957.63824
74 3384.03824 607.53824
75 628.03824 3384.03824
76 1856.53824 628.03824
77 2877.83824 1856.53824
78 1673.63824 2877.83824
79 -231.16176 1673.63824
80 1847.03824 -231.16176
81 4193.43824 1847.03824
82 2697.73824 4193.43824
83 -89.16176 2697.73824
84 3244.63824 -89.16176
85 3456.53824 3244.63824
86 3722.13824 3456.53824
87 5276.63824 3722.13824
88 3523.63824 5276.63824
89 4959.23824 3523.63824
90 5218.23824 4959.23824
91 595.63824 5218.23824
> 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/7wvta1229016224.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/8uc2u1229016224.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/9thre1229016224.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/10c07a1229016224.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/11v7xa1229016224.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/12rmi61229016224.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/13l9s11229016225.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/14hxxq1229016225.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/15nnf41229016225.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/16a4xi1229016225.tab")
+ }
>
> system("convert tmp/19ulz1229016224.ps tmp/19ulz1229016224.png")
> system("convert tmp/2fea61229016224.ps tmp/2fea61229016224.png")
> system("convert tmp/39m7w1229016224.ps tmp/39m7w1229016224.png")
> system("convert tmp/4xime1229016224.ps tmp/4xime1229016224.png")
> system("convert tmp/5v1kj1229016224.ps tmp/5v1kj1229016224.png")
> system("convert tmp/6vdrs1229016224.ps tmp/6vdrs1229016224.png")
> system("convert tmp/7wvta1229016224.ps tmp/7wvta1229016224.png")
> system("convert tmp/8uc2u1229016224.ps tmp/8uc2u1229016224.png")
> system("convert tmp/9thre1229016224.ps tmp/9thre1229016224.png")
> system("convert tmp/10c07a1229016224.ps tmp/10c07a1229016224.png")
>
>
> proc.time()
user system elapsed
2.826 1.651 3.270