R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(100.00,0,100.42,0,100.50,0,101.14,0,101.98,0,102.31,0,103.27,0,103.80,0,103.46,0,105.06,0,106.08,0,106.74,0,107.35,0,108.96,0,109.85,0,109.81,0,109.99,0,111.60,0,112.74,0,112.78,0,113.66,0,115.37,0,116.26,0,116.24,0,116.73,0,118.76,0,119.78,0,120.23,0,121.48,0,124.07,0,125.82,0,126.92,0,128.48,0,131.44,0,133.51,0,134.58,0,136.68,0,140.10,0,142.45,0,143.91,0,146.19,0,149.84,0,152.31,0,153.62,0,155.79,0,159.89,0,163.21,0,165.32,0,167.68,0,171.79,0,175.38,0,177.81,0,181.09,0,186.48,0,191.07,0,194.23,0,197.82,0,204.41,0,209.26,0,212.24,0,214.88,0,218.87,0,219.86,0,219.75,0,220.89,0,224.02,0,222.27,0,217.27,1,213.23,1,212.44,1,207.87,1,199.46,1,198.19,1,199.77,1,200.10,1,195.76,1,191.27,1,195.79,1,192.7,1),dim=c(2,79),dimnames=list(c('woningprijsindex_us','Dummy_'),1:79))
> y <- array(NA,dim=c(2,79),dimnames=list(c('woningprijsindex_us','Dummy_'),1:79))
> 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'
> #'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
> 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
woningprijsindex_us Dummy_ t
1 100.00 0 1
2 100.42 0 2
3 100.50 0 3
4 101.14 0 4
5 101.98 0 5
6 102.31 0 6
7 103.27 0 7
8 103.80 0 8
9 103.46 0 9
10 105.06 0 10
11 106.08 0 11
12 106.74 0 12
13 107.35 0 13
14 108.96 0 14
15 109.85 0 15
16 109.81 0 16
17 109.99 0 17
18 111.60 0 18
19 112.74 0 19
20 112.78 0 20
21 113.66 0 21
22 115.37 0 22
23 116.26 0 23
24 116.24 0 24
25 116.73 0 25
26 118.76 0 26
27 119.78 0 27
28 120.23 0 28
29 121.48 0 29
30 124.07 0 30
31 125.82 0 31
32 126.92 0 32
33 128.48 0 33
34 131.44 0 34
35 133.51 0 35
36 134.58 0 36
37 136.68 0 37
38 140.10 0 38
39 142.45 0 39
40 143.91 0 40
41 146.19 0 41
42 149.84 0 42
43 152.31 0 43
44 153.62 0 44
45 155.79 0 45
46 159.89 0 46
47 163.21 0 47
48 165.32 0 48
49 167.68 0 49
50 171.79 0 50
51 175.38 0 51
52 177.81 0 52
53 181.09 0 53
54 186.48 0 54
55 191.07 0 55
56 194.23 0 56
57 197.82 0 57
58 204.41 0 58
59 209.26 0 59
60 212.24 0 60
61 214.88 0 61
62 218.87 0 62
63 219.86 0 63
64 219.75 0 64
65 220.89 0 65
66 224.02 0 66
67 222.27 0 67
68 217.27 1 68
69 213.23 1 69
70 212.44 1 70
71 207.87 1 71
72 199.46 1 72
73 198.19 1 73
74 199.77 1 74
75 200.10 1 75
76 195.76 1 76
77 191.27 1 77
78 195.79 1 78
79 192.70 1 79
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Dummy_ t
78.444 -19.860 1.951
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-20.018 -10.610 -2.568 10.964 26.013
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 78.444 2.975 26.371 < 2e-16 ***
Dummy_ -19.860 4.828 -4.113 9.8e-05 ***
t 1.951 0.076 25.673 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 12.06 on 76 degrees of freedom
Multiple R-squared: 0.9212, Adjusted R-squared: 0.9191
F-statistic: 444 on 2 and 76 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,] 1.646698e-05 3.293395e-05 9.999835e-01
[2,] 9.995052e-07 1.999010e-06 9.999990e-01
[3,] 3.292642e-08 6.585284e-08 1.000000e+00
[4,] 5.352165e-09 1.070433e-08 1.000000e+00
[5,] 4.978428e-10 9.956855e-10 1.000000e+00
[6,] 1.025228e-10 2.050456e-10 1.000000e+00
[7,] 1.080259e-11 2.160517e-11 1.000000e+00
[8,] 7.875849e-13 1.575170e-12 1.000000e+00
[9,] 6.185067e-13 1.237013e-12 1.000000e+00
[10,] 2.226364e-13 4.452727e-13 1.000000e+00
[11,] 1.648735e-14 3.297471e-14 1.000000e+00
[12,] 1.355135e-15 2.710270e-15 1.000000e+00
[13,] 1.553325e-16 3.106650e-16 1.000000e+00
[14,] 3.276464e-17 6.552927e-17 1.000000e+00
[15,] 2.435095e-18 4.870190e-18 1.000000e+00
[16,] 1.866196e-19 3.732391e-19 1.000000e+00
[17,] 7.720676e-20 1.544135e-19 1.000000e+00
[18,] 2.540048e-20 5.080096e-20 1.000000e+00
[19,] 2.012892e-21 4.025784e-21 1.000000e+00
[20,] 1.466497e-22 2.932994e-22 1.000000e+00
[21,] 5.309816e-23 1.061963e-22 1.000000e+00
[22,] 2.126040e-23 4.252079e-23 1.000000e+00
[23,] 2.795488e-24 5.590977e-24 1.000000e+00
[24,] 7.020188e-25 1.404038e-24 1.000000e+00
[25,] 2.073906e-23 4.147812e-23 1.000000e+00
[26,] 6.638579e-22 1.327716e-21 1.000000e+00
[27,] 3.566164e-21 7.132328e-21 1.000000e+00
[28,] 1.790247e-20 3.580495e-20 1.000000e+00
[29,] 9.415218e-19 1.883044e-18 1.000000e+00
[30,] 2.696870e-17 5.393741e-17 1.000000e+00
[31,] 1.396416e-16 2.792833e-16 1.000000e+00
[32,] 8.153583e-16 1.630717e-15 1.000000e+00
[33,] 1.700527e-14 3.401053e-14 1.000000e+00
[34,] 2.325785e-13 4.651570e-13 1.000000e+00
[35,] 1.209820e-12 2.419639e-12 1.000000e+00
[36,] 5.766566e-12 1.153313e-11 1.000000e+00
[37,] 5.232623e-11 1.046525e-10 1.000000e+00
[38,] 3.521559e-10 7.043119e-10 1.000000e+00
[39,] 1.217720e-09 2.435439e-09 1.000000e+00
[40,] 3.911171e-09 7.822341e-09 1.000000e+00
[41,] 2.115563e-08 4.231126e-08 1.000000e+00
[42,] 1.236719e-07 2.473439e-07 9.999999e-01
[43,] 5.679938e-07 1.135988e-06 9.999994e-01
[44,] 2.537842e-06 5.075685e-06 9.999975e-01
[45,] 1.454341e-05 2.908682e-05 9.999855e-01
[46,] 8.865320e-05 1.773064e-04 9.999113e-01
[47,] 5.299012e-04 1.059802e-03 9.994701e-01
[48,] 3.458978e-03 6.917956e-03 9.965410e-01
[49,] 2.130458e-02 4.260916e-02 9.786954e-01
[50,] 9.872206e-02 1.974441e-01 9.012779e-01
[51,] 3.259099e-01 6.518197e-01 6.740901e-01
[52,] 7.038351e-01 5.923298e-01 2.961649e-01
[53,] 9.257740e-01 1.484520e-01 7.422602e-02
[54,] 9.868766e-01 2.624674e-02 1.312337e-02
[55,] 9.982713e-01 3.457312e-03 1.728656e-03
[56,] 9.997990e-01 4.020318e-04 2.010159e-04
[57,] 9.999305e-01 1.389310e-04 6.946549e-05
[58,] 9.999594e-01 8.128146e-05 4.064073e-05
[59,] 9.999727e-01 5.458739e-05 2.729370e-05
[60,] 9.999661e-01 6.785522e-05 3.392761e-05
[61,] 9.998951e-01 2.098407e-04 1.049203e-04
[62,] 9.996553e-01 6.894301e-04 3.447151e-04
[63,] 9.992757e-01 1.448698e-03 7.243491e-04
[64,] 9.981740e-01 3.651974e-03 1.825987e-03
[65,] 9.983088e-01 3.382303e-03 1.691152e-03
[66,] 9.985097e-01 2.980589e-03 1.490295e-03
[67,] 9.937874e-01 1.242528e-02 6.212638e-03
[68,] 9.792267e-01 4.154656e-02 2.077328e-02
> postscript(file="/var/www/rcomp/tmp/1f1gj1292593918.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/2f1gj1292593918.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/38sy41292593918.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/48sy41292593918.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/58sy41292593918.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 = 79
Frequency = 1
1 2 3 4 5 6
19.6047495 18.0736874 16.2026254 14.8915633 13.7805013 12.1594392
7 8 9 10 11 12
11.1683771 9.7473151 7.4562530 7.1051910 6.1741289 4.8830668
13 14 15 16 17 18
3.5420048 3.2009427 2.1398807 0.1488186 -1.6222435 -1.9633055
19 20 21 22 23 24
-2.7743676 -4.6854297 -5.7564917 -5.9975538 -7.0586158 -9.0296779
25 26 27 28 29 30
-10.4907400 -10.4118020 -11.3428641 -12.8439261 -13.5449882 -12.9060503
31 32 33 34 35 36
-13.1071123 -13.9581744 -14.3492364 -13.3402985 -13.2213606 -14.1024226
37 38 39 40 41 42
-13.9534847 -12.4845468 -12.0856088 -12.5766709 -12.2477329 -10.5487950
43 44 45 46 47 48
-10.0298571 -10.6709191 -10.4519812 -8.3030432 -6.9341053 -6.7751674
49 50 51 52 53 54
-6.3662294 -4.2072915 -2.5683535 -2.0894156 -0.7604777 2.6784603
55 56 57 58 59 60
5.3173982 6.5263361 8.1652741 12.8042120 15.7031500 16.7320879
61 62 63 64 65 66
17.4210258 19.4599638 18.4989017 16.4378397 15.6267776 16.8057155
67 68 69 70 71 72
13.1046535 26.0133413 20.0222793 17.2812172 10.7601552 0.3990931
73 74 75 76 77 78
-2.8219690 -3.1930310 -4.8140931 -11.1051552 -17.5462172 -14.9772793
79
-20.0183413
> postscript(file="/var/www/rcomp/tmp/61jf71292593918.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 = 79
Frequency = 1
lag(myerror, k = 1) myerror
0 19.6047495 NA
1 18.0736874 19.6047495
2 16.2026254 18.0736874
3 14.8915633 16.2026254
4 13.7805013 14.8915633
5 12.1594392 13.7805013
6 11.1683771 12.1594392
7 9.7473151 11.1683771
8 7.4562530 9.7473151
9 7.1051910 7.4562530
10 6.1741289 7.1051910
11 4.8830668 6.1741289
12 3.5420048 4.8830668
13 3.2009427 3.5420048
14 2.1398807 3.2009427
15 0.1488186 2.1398807
16 -1.6222435 0.1488186
17 -1.9633055 -1.6222435
18 -2.7743676 -1.9633055
19 -4.6854297 -2.7743676
20 -5.7564917 -4.6854297
21 -5.9975538 -5.7564917
22 -7.0586158 -5.9975538
23 -9.0296779 -7.0586158
24 -10.4907400 -9.0296779
25 -10.4118020 -10.4907400
26 -11.3428641 -10.4118020
27 -12.8439261 -11.3428641
28 -13.5449882 -12.8439261
29 -12.9060503 -13.5449882
30 -13.1071123 -12.9060503
31 -13.9581744 -13.1071123
32 -14.3492364 -13.9581744
33 -13.3402985 -14.3492364
34 -13.2213606 -13.3402985
35 -14.1024226 -13.2213606
36 -13.9534847 -14.1024226
37 -12.4845468 -13.9534847
38 -12.0856088 -12.4845468
39 -12.5766709 -12.0856088
40 -12.2477329 -12.5766709
41 -10.5487950 -12.2477329
42 -10.0298571 -10.5487950
43 -10.6709191 -10.0298571
44 -10.4519812 -10.6709191
45 -8.3030432 -10.4519812
46 -6.9341053 -8.3030432
47 -6.7751674 -6.9341053
48 -6.3662294 -6.7751674
49 -4.2072915 -6.3662294
50 -2.5683535 -4.2072915
51 -2.0894156 -2.5683535
52 -0.7604777 -2.0894156
53 2.6784603 -0.7604777
54 5.3173982 2.6784603
55 6.5263361 5.3173982
56 8.1652741 6.5263361
57 12.8042120 8.1652741
58 15.7031500 12.8042120
59 16.7320879 15.7031500
60 17.4210258 16.7320879
61 19.4599638 17.4210258
62 18.4989017 19.4599638
63 16.4378397 18.4989017
64 15.6267776 16.4378397
65 16.8057155 15.6267776
66 13.1046535 16.8057155
67 26.0133413 13.1046535
68 20.0222793 26.0133413
69 17.2812172 20.0222793
70 10.7601552 17.2812172
71 0.3990931 10.7601552
72 -2.8219690 0.3990931
73 -3.1930310 -2.8219690
74 -4.8140931 -3.1930310
75 -11.1051552 -4.8140931
76 -17.5462172 -11.1051552
77 -14.9772793 -17.5462172
78 -20.0183413 -14.9772793
79 NA -20.0183413
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 18.0736874 19.6047495
[2,] 16.2026254 18.0736874
[3,] 14.8915633 16.2026254
[4,] 13.7805013 14.8915633
[5,] 12.1594392 13.7805013
[6,] 11.1683771 12.1594392
[7,] 9.7473151 11.1683771
[8,] 7.4562530 9.7473151
[9,] 7.1051910 7.4562530
[10,] 6.1741289 7.1051910
[11,] 4.8830668 6.1741289
[12,] 3.5420048 4.8830668
[13,] 3.2009427 3.5420048
[14,] 2.1398807 3.2009427
[15,] 0.1488186 2.1398807
[16,] -1.6222435 0.1488186
[17,] -1.9633055 -1.6222435
[18,] -2.7743676 -1.9633055
[19,] -4.6854297 -2.7743676
[20,] -5.7564917 -4.6854297
[21,] -5.9975538 -5.7564917
[22,] -7.0586158 -5.9975538
[23,] -9.0296779 -7.0586158
[24,] -10.4907400 -9.0296779
[25,] -10.4118020 -10.4907400
[26,] -11.3428641 -10.4118020
[27,] -12.8439261 -11.3428641
[28,] -13.5449882 -12.8439261
[29,] -12.9060503 -13.5449882
[30,] -13.1071123 -12.9060503
[31,] -13.9581744 -13.1071123
[32,] -14.3492364 -13.9581744
[33,] -13.3402985 -14.3492364
[34,] -13.2213606 -13.3402985
[35,] -14.1024226 -13.2213606
[36,] -13.9534847 -14.1024226
[37,] -12.4845468 -13.9534847
[38,] -12.0856088 -12.4845468
[39,] -12.5766709 -12.0856088
[40,] -12.2477329 -12.5766709
[41,] -10.5487950 -12.2477329
[42,] -10.0298571 -10.5487950
[43,] -10.6709191 -10.0298571
[44,] -10.4519812 -10.6709191
[45,] -8.3030432 -10.4519812
[46,] -6.9341053 -8.3030432
[47,] -6.7751674 -6.9341053
[48,] -6.3662294 -6.7751674
[49,] -4.2072915 -6.3662294
[50,] -2.5683535 -4.2072915
[51,] -2.0894156 -2.5683535
[52,] -0.7604777 -2.0894156
[53,] 2.6784603 -0.7604777
[54,] 5.3173982 2.6784603
[55,] 6.5263361 5.3173982
[56,] 8.1652741 6.5263361
[57,] 12.8042120 8.1652741
[58,] 15.7031500 12.8042120
[59,] 16.7320879 15.7031500
[60,] 17.4210258 16.7320879
[61,] 19.4599638 17.4210258
[62,] 18.4989017 19.4599638
[63,] 16.4378397 18.4989017
[64,] 15.6267776 16.4378397
[65,] 16.8057155 15.6267776
[66,] 13.1046535 16.8057155
[67,] 26.0133413 13.1046535
[68,] 20.0222793 26.0133413
[69,] 17.2812172 20.0222793
[70,] 10.7601552 17.2812172
[71,] 0.3990931 10.7601552
[72,] -2.8219690 0.3990931
[73,] -3.1930310 -2.8219690
[74,] -4.8140931 -3.1930310
[75,] -11.1051552 -4.8140931
[76,] -17.5462172 -11.1051552
[77,] -14.9772793 -17.5462172
[78,] -20.0183413 -14.9772793
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 18.0736874 19.6047495
2 16.2026254 18.0736874
3 14.8915633 16.2026254
4 13.7805013 14.8915633
5 12.1594392 13.7805013
6 11.1683771 12.1594392
7 9.7473151 11.1683771
8 7.4562530 9.7473151
9 7.1051910 7.4562530
10 6.1741289 7.1051910
11 4.8830668 6.1741289
12 3.5420048 4.8830668
13 3.2009427 3.5420048
14 2.1398807 3.2009427
15 0.1488186 2.1398807
16 -1.6222435 0.1488186
17 -1.9633055 -1.6222435
18 -2.7743676 -1.9633055
19 -4.6854297 -2.7743676
20 -5.7564917 -4.6854297
21 -5.9975538 -5.7564917
22 -7.0586158 -5.9975538
23 -9.0296779 -7.0586158
24 -10.4907400 -9.0296779
25 -10.4118020 -10.4907400
26 -11.3428641 -10.4118020
27 -12.8439261 -11.3428641
28 -13.5449882 -12.8439261
29 -12.9060503 -13.5449882
30 -13.1071123 -12.9060503
31 -13.9581744 -13.1071123
32 -14.3492364 -13.9581744
33 -13.3402985 -14.3492364
34 -13.2213606 -13.3402985
35 -14.1024226 -13.2213606
36 -13.9534847 -14.1024226
37 -12.4845468 -13.9534847
38 -12.0856088 -12.4845468
39 -12.5766709 -12.0856088
40 -12.2477329 -12.5766709
41 -10.5487950 -12.2477329
42 -10.0298571 -10.5487950
43 -10.6709191 -10.0298571
44 -10.4519812 -10.6709191
45 -8.3030432 -10.4519812
46 -6.9341053 -8.3030432
47 -6.7751674 -6.9341053
48 -6.3662294 -6.7751674
49 -4.2072915 -6.3662294
50 -2.5683535 -4.2072915
51 -2.0894156 -2.5683535
52 -0.7604777 -2.0894156
53 2.6784603 -0.7604777
54 5.3173982 2.6784603
55 6.5263361 5.3173982
56 8.1652741 6.5263361
57 12.8042120 8.1652741
58 15.7031500 12.8042120
59 16.7320879 15.7031500
60 17.4210258 16.7320879
61 19.4599638 17.4210258
62 18.4989017 19.4599638
63 16.4378397 18.4989017
64 15.6267776 16.4378397
65 16.8057155 15.6267776
66 13.1046535 16.8057155
67 26.0133413 13.1046535
68 20.0222793 26.0133413
69 17.2812172 20.0222793
70 10.7601552 17.2812172
71 0.3990931 10.7601552
72 -2.8219690 0.3990931
73 -3.1930310 -2.8219690
74 -4.8140931 -3.1930310
75 -11.1051552 -4.8140931
76 -17.5462172 -11.1051552
77 -14.9772793 -17.5462172
78 -20.0183413 -14.9772793
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/7usws1292593918.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/8usws1292593918.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/rcomp/tmp/9usws1292593918.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/rcomp/tmp/1042vv1292593918.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/1182u01292593918.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/12blbp1292593918.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/1304701292593918.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/rcomp/tmp/14sd7l1292593918.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/rcomp/tmp/15ww591292593918.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/rcomp/tmp/16snli1292593918.tab")
+ }
>
> try(system("convert tmp/1f1gj1292593918.ps tmp/1f1gj1292593918.png",intern=TRUE))
character(0)
> try(system("convert tmp/2f1gj1292593918.ps tmp/2f1gj1292593918.png",intern=TRUE))
character(0)
> try(system("convert tmp/38sy41292593918.ps tmp/38sy41292593918.png",intern=TRUE))
character(0)
> try(system("convert tmp/48sy41292593918.ps tmp/48sy41292593918.png",intern=TRUE))
character(0)
> try(system("convert tmp/58sy41292593918.ps tmp/58sy41292593918.png",intern=TRUE))
character(0)
> try(system("convert tmp/61jf71292593918.ps tmp/61jf71292593918.png",intern=TRUE))
character(0)
> try(system("convert tmp/7usws1292593918.ps tmp/7usws1292593918.png",intern=TRUE))
character(0)
> try(system("convert tmp/8usws1292593918.ps tmp/8usws1292593918.png",intern=TRUE))
character(0)
> try(system("convert tmp/9usws1292593918.ps tmp/9usws1292593918.png",intern=TRUE))
character(0)
> try(system("convert tmp/1042vv1292593918.ps tmp/1042vv1292593918.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.260 1.720 4.971