R version 2.13.0 (2011-04-13)
Copyright (C) 2011 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(465000
+ ,1520
+ ,510
+ ,3
+ ,1979
+ ,530000
+ ,2700
+ ,345
+ ,4
+ ,1977
+ ,389500
+ ,3571
+ ,150
+ ,4
+ ,1969
+ ,305000
+ ,854
+ ,260
+ ,3
+ ,2011
+ ,620000
+ ,1560
+ ,458
+ ,5
+ ,1981
+ ,750000
+ ,3017
+ ,400
+ ,4
+ ,1972
+ ,389000
+ ,436
+ ,201
+ ,4
+ ,2011
+ ,387000
+ ,1098
+ ,233
+ ,4
+ ,1966
+ ,312000
+ ,625
+ ,160
+ ,2
+ ,1953
+ ,375000
+ ,700
+ ,140
+ ,3
+ ,1985
+ ,385000
+ ,639
+ ,155
+ ,2
+ ,1990
+ ,395000
+ ,1246
+ ,300
+ ,5
+ ,1973
+ ,398000
+ ,600
+ ,220
+ ,3
+ ,1997
+ ,449000
+ ,1000
+ ,272
+ ,4
+ ,1982
+ ,451245
+ ,1047
+ ,280
+ ,4
+ ,2011
+ ,511862
+ ,1414
+ ,219
+ ,3
+ ,2011
+ ,324000
+ ,674
+ ,160
+ ,3
+ ,1973
+ ,772000
+ ,6500
+ ,190
+ ,5
+ ,1971
+ ,617000
+ ,3321
+ ,192
+ ,3
+ ,1980
+ ,595000
+ ,2000
+ ,320
+ ,5
+ ,1997
+ ,475000
+ ,1945
+ ,241
+ ,4
+ ,1974
+ ,985000
+ ,7620
+ ,500
+ ,4
+ ,1977
+ ,439000
+ ,1001
+ ,190
+ ,3
+ ,1991
+ ,479000
+ ,1699
+ ,261
+ ,3
+ ,1987
+ ,657160
+ ,1961
+ ,206
+ ,3
+ ,2012
+ ,299000
+ ,1248
+ ,127
+ ,1
+ ,1970
+ ,419000
+ ,1294
+ ,225
+ ,4
+ ,1989
+ ,449000
+ ,1267
+ ,185
+ ,4
+ ,1970
+ ,327000
+ ,998
+ ,215
+ ,4
+ ,1990
+ ,1695000
+ ,5462
+ ,730
+ ,4
+ ,1998
+ ,489000
+ ,1883
+ ,223
+ ,3
+ ,1987
+ ,449000
+ ,1000
+ ,256
+ ,4
+ ,1994
+ ,470000
+ ,663
+ ,281
+ ,4
+ ,2011
+ ,537000
+ ,2240
+ ,298
+ ,3
+ ,1976
+ ,685000
+ ,2580
+ ,362
+ ,4
+ ,2001
+ ,399000
+ ,2755
+ ,250
+ ,3
+ ,1980
+ ,299500
+ ,773
+ ,188
+ ,4
+ ,1968
+ ,598000
+ ,1465
+ ,500
+ ,4
+ ,1982
+ ,547000
+ ,2025
+ ,270
+ ,4
+ ,1965
+ ,750000
+ ,2160
+ ,300
+ ,5
+ ,1961
+ ,320000
+ ,983
+ ,130
+ ,2
+ ,1989
+ ,373000
+ ,351
+ ,200
+ ,4
+ ,2011
+ ,825000
+ ,712
+ ,270
+ ,3
+ ,2003
+ ,389000
+ ,1120
+ ,224
+ ,4
+ ,1975
+ ,474000
+ ,2619
+ ,290
+ ,4
+ ,1967
+ ,325000
+ ,1193
+ ,214
+ ,3
+ ,1964
+ ,795000
+ ,1500
+ ,450
+ ,4
+ ,2011
+ ,590000
+ ,8560
+ ,330
+ ,3
+ ,1950
+ ,608000
+ ,2236
+ ,190
+ ,3
+ ,1993
+ ,1300000
+ ,3390
+ ,462
+ ,3
+ ,1993
+ ,1325000
+ ,2935
+ ,473
+ ,4
+ ,2004
+ ,1680000
+ ,3700
+ ,528
+ ,3
+ ,2008
+ ,895000
+ ,3290
+ ,470
+ ,6
+ ,1973
+ ,235000
+ ,1115
+ ,94
+ ,2
+ ,1938
+ ,330000
+ ,1200
+ ,100
+ ,3
+ ,1970
+ ,489000
+ ,2160
+ ,166
+ ,4
+ ,1977
+ ,499000
+ ,2605
+ ,334
+ ,5
+ ,1963
+ ,535000
+ ,2229
+ ,230
+ ,4
+ ,1974
+ ,645000
+ ,2267
+ ,303
+ ,4
+ ,1980
+ ,699000
+ ,5027
+ ,315
+ ,5
+ ,1976)
+ ,dim=c(5
+ ,60)
+ ,dimnames=list(c('Prijs'
+ ,'Domeinopp.'
+ ,'Huisopp.'
+ ,'Slaapkamers'
+ ,'bouwjaar')
+ ,1:60))
> y <- array(NA,dim=c(5,60),dimnames=list(c('Prijs','Domeinopp.','Huisopp.','Slaapkamers','bouwjaar'),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 = 'Do not include Seasonal Dummies'
> par1 = '1'
> 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
Prijs Domeinopp. Huisopp. Slaapkamers bouwjaar
1 465000 1520 510 3 1979
2 530000 2700 345 4 1977
3 389500 3571 150 4 1969
4 305000 854 260 3 2011
5 620000 1560 458 5 1981
6 750000 3017 400 4 1972
7 389000 436 201 4 2011
8 387000 1098 233 4 1966
9 312000 625 160 2 1953
10 375000 700 140 3 1985
11 385000 639 155 2 1990
12 395000 1246 300 5 1973
13 398000 600 220 3 1997
14 449000 1000 272 4 1982
15 451245 1047 280 4 2011
16 511862 1414 219 3 2011
17 324000 674 160 3 1973
18 772000 6500 190 5 1971
19 617000 3321 192 3 1980
20 595000 2000 320 5 1997
21 475000 1945 241 4 1974
22 985000 7620 500 4 1977
23 439000 1001 190 3 1991
24 479000 1699 261 3 1987
25 657160 1961 206 3 2012
26 299000 1248 127 1 1970
27 419000 1294 225 4 1989
28 449000 1267 185 4 1970
29 327000 998 215 4 1990
30 1695000 5462 730 4 1998
31 489000 1883 223 3 1987
32 449000 1000 256 4 1994
33 470000 663 281 4 2011
34 537000 2240 298 3 1976
35 685000 2580 362 4 2001
36 399000 2755 250 3 1980
37 299500 773 188 4 1968
38 598000 1465 500 4 1982
39 547000 2025 270 4 1965
40 750000 2160 300 5 1961
41 320000 983 130 2 1989
42 373000 351 200 4 2011
43 825000 712 270 3 2003
44 389000 1120 224 4 1975
45 474000 2619 290 4 1967
46 325000 1193 214 3 1964
47 795000 1500 450 4 2011
48 590000 8560 330 3 1950
49 608000 2236 190 3 1993
50 1300000 3390 462 3 1993
51 1325000 2935 473 4 2004
52 1680000 3700 528 3 2008
53 895000 3290 470 6 1973
54 235000 1115 94 2 1938
55 330000 1200 100 3 1970
56 489000 2160 166 4 1977
57 499000 2605 334 5 1963
58 535000 2229 230 4 1974
59 645000 2267 303 4 1980
60 699000 5027 315 5 1976
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Domeinopp. Huisopp. Slaapkamers bouwjaar
-8.236e+06 6.474e+01 1.460e+03 -3.361e+04 4.228e+03
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-408464 -69867 -22530 83926 516507
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8.236e+06 2.743e+06 -3.002 0.004023 **
Domeinopp. 6.474e+01 1.561e+01 4.147 0.000117 ***
Huisopp. 1.460e+03 2.194e+02 6.656 1.38e-08 ***
Slaapkamers -3.361e+04 2.594e+04 -1.296 0.200508
bouwjaar 4.228e+03 1.384e+03 3.054 0.003476 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 164200 on 55 degrees of freedom
Multiple R-squared: 0.7229, Adjusted R-squared: 0.7027
F-statistic: 35.87 on 4 and 55 DF, p-value: 9.826e-15
> 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,] 2.767506e-01 5.535012e-01 0.72324939
[2,] 2.238539e-01 4.477078e-01 0.77614612
[3,] 1.556153e-01 3.112306e-01 0.84438472
[4,] 1.283660e-01 2.567320e-01 0.87163399
[5,] 9.255119e-02 1.851024e-01 0.90744881
[6,] 5.268828e-02 1.053766e-01 0.94731172
[7,] 2.757954e-02 5.515908e-02 0.97242046
[8,] 1.524323e-02 3.048646e-02 0.98475677
[9,] 1.121668e-02 2.243336e-02 0.98878332
[10,] 5.237243e-03 1.047449e-02 0.99476276
[11,] 3.998872e-03 7.997744e-03 0.99600113
[12,] 2.911882e-03 5.823764e-03 0.99708812
[13,] 1.625266e-03 3.250532e-03 0.99837473
[14,] 7.084238e-04 1.416848e-03 0.99929158
[15,] 3.822129e-04 7.644257e-04 0.99961779
[16,] 1.825564e-04 3.651129e-04 0.99981744
[17,] 8.170505e-05 1.634101e-04 0.99991829
[18,] 1.533461e-04 3.066921e-04 0.99984665
[19,] 7.617099e-05 1.523420e-04 0.99992383
[20,] 3.186844e-05 6.373688e-05 0.99996813
[21,] 1.993630e-05 3.987260e-05 0.99998006
[22,] 1.558032e-05 3.116065e-05 0.99998442
[23,] 2.581531e-02 5.163062e-02 0.97418469
[24,] 1.601337e-02 3.202675e-02 0.98398663
[25,] 1.009583e-02 2.019167e-02 0.98990417
[26,] 8.251934e-03 1.650387e-02 0.99174807
[27,] 5.395618e-03 1.079124e-02 0.99460438
[28,] 4.216222e-03 8.432444e-03 0.99578378
[29,] 6.291649e-03 1.258330e-02 0.99370835
[30,] 3.496795e-03 6.993590e-03 0.99650321
[31,] 3.389011e-02 6.778023e-02 0.96610989
[32,] 2.335865e-02 4.671729e-02 0.97664135
[33,] 6.946763e-02 1.389353e-01 0.93053237
[34,] 6.245487e-02 1.249097e-01 0.93754513
[35,] 7.165691e-02 1.433138e-01 0.92834309
[36,] 1.204810e-01 2.409620e-01 0.87951898
[37,] 8.504472e-02 1.700894e-01 0.91495528
[38,] 6.521771e-02 1.304354e-01 0.93478229
[39,] 6.303582e-02 1.260716e-01 0.93696418
[40,] 8.729185e-01 2.541631e-01 0.12708155
[41,] 9.294731e-01 1.410538e-01 0.07052691
[42,] 9.324360e-01 1.351280e-01 0.06756401
[43,] 9.372167e-01 1.255665e-01 0.06278327
[44,] 8.939459e-01 2.121081e-01 0.10605406
[45,] 8.974807e-01 2.050387e-01 0.10251933
> postscript(file="/var/wessaorg/rcomp/tmp/1nwfn1324475533.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/wessaorg/rcomp/tmp/28q8c1324475533.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/wessaorg/rcomp/tmp/399cq1324475533.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/wessaorg/rcomp/tmp/4agti1324475533.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/wessaorg/rcomp/tmp/5pv3k1324475533.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 = 60
Frequency = 1
1 2 3 4 5 6
-408463.718 -136858.130 -15181.566 -295596.702 -121363.082 3449.126
7 8 9 10 11 12
-64776.050 33905.086 83870.631 69527.901 6825.466 -61498.490
13 14 15 16 17 18
-68551.553 -22348.013 -157442.345 -65120.179 41744.328 144442.219
19 20 21 22 23 24
87057.579 -40990.047 21564.890 -226704.177 15662.846 -76286.701
25 26 27 28 29 30
59520.212 -26760.725 -32348.488 138141.472 -94811.896 198365.702
31 32 33 34 35 36
-22710.995 -49722.035 -115287.849 -60828.694 -100387.377 -178991.931
37 38 39 40 41 42
24698.039 -236377.416 84092.884 285067.493 -39711.547 -73813.051
43 44 45 46 47 48
252819.023 9569.712 -65022.190 -31653.038 -91248.216 -353773.621
49 50 51 52 53 54
96254.228 316370.944 341863.712 516507.090 91549.184 134943.450
55 56 57 58 59 60
113988.102 118476.826 -52844.317 79241.286 54817.906 -36863.199
> postscript(file="/var/wessaorg/rcomp/tmp/6jxgz1324475533.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -408463.718 NA
1 -136858.130 -408463.718
2 -15181.566 -136858.130
3 -295596.702 -15181.566
4 -121363.082 -295596.702
5 3449.126 -121363.082
6 -64776.050 3449.126
7 33905.086 -64776.050
8 83870.631 33905.086
9 69527.901 83870.631
10 6825.466 69527.901
11 -61498.490 6825.466
12 -68551.553 -61498.490
13 -22348.013 -68551.553
14 -157442.345 -22348.013
15 -65120.179 -157442.345
16 41744.328 -65120.179
17 144442.219 41744.328
18 87057.579 144442.219
19 -40990.047 87057.579
20 21564.890 -40990.047
21 -226704.177 21564.890
22 15662.846 -226704.177
23 -76286.701 15662.846
24 59520.212 -76286.701
25 -26760.725 59520.212
26 -32348.488 -26760.725
27 138141.472 -32348.488
28 -94811.896 138141.472
29 198365.702 -94811.896
30 -22710.995 198365.702
31 -49722.035 -22710.995
32 -115287.849 -49722.035
33 -60828.694 -115287.849
34 -100387.377 -60828.694
35 -178991.931 -100387.377
36 24698.039 -178991.931
37 -236377.416 24698.039
38 84092.884 -236377.416
39 285067.493 84092.884
40 -39711.547 285067.493
41 -73813.051 -39711.547
42 252819.023 -73813.051
43 9569.712 252819.023
44 -65022.190 9569.712
45 -31653.038 -65022.190
46 -91248.216 -31653.038
47 -353773.621 -91248.216
48 96254.228 -353773.621
49 316370.944 96254.228
50 341863.712 316370.944
51 516507.090 341863.712
52 91549.184 516507.090
53 134943.450 91549.184
54 113988.102 134943.450
55 118476.826 113988.102
56 -52844.317 118476.826
57 79241.286 -52844.317
58 54817.906 79241.286
59 -36863.199 54817.906
60 NA -36863.199
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -136858.130 -408463.718
[2,] -15181.566 -136858.130
[3,] -295596.702 -15181.566
[4,] -121363.082 -295596.702
[5,] 3449.126 -121363.082
[6,] -64776.050 3449.126
[7,] 33905.086 -64776.050
[8,] 83870.631 33905.086
[9,] 69527.901 83870.631
[10,] 6825.466 69527.901
[11,] -61498.490 6825.466
[12,] -68551.553 -61498.490
[13,] -22348.013 -68551.553
[14,] -157442.345 -22348.013
[15,] -65120.179 -157442.345
[16,] 41744.328 -65120.179
[17,] 144442.219 41744.328
[18,] 87057.579 144442.219
[19,] -40990.047 87057.579
[20,] 21564.890 -40990.047
[21,] -226704.177 21564.890
[22,] 15662.846 -226704.177
[23,] -76286.701 15662.846
[24,] 59520.212 -76286.701
[25,] -26760.725 59520.212
[26,] -32348.488 -26760.725
[27,] 138141.472 -32348.488
[28,] -94811.896 138141.472
[29,] 198365.702 -94811.896
[30,] -22710.995 198365.702
[31,] -49722.035 -22710.995
[32,] -115287.849 -49722.035
[33,] -60828.694 -115287.849
[34,] -100387.377 -60828.694
[35,] -178991.931 -100387.377
[36,] 24698.039 -178991.931
[37,] -236377.416 24698.039
[38,] 84092.884 -236377.416
[39,] 285067.493 84092.884
[40,] -39711.547 285067.493
[41,] -73813.051 -39711.547
[42,] 252819.023 -73813.051
[43,] 9569.712 252819.023
[44,] -65022.190 9569.712
[45,] -31653.038 -65022.190
[46,] -91248.216 -31653.038
[47,] -353773.621 -91248.216
[48,] 96254.228 -353773.621
[49,] 316370.944 96254.228
[50,] 341863.712 316370.944
[51,] 516507.090 341863.712
[52,] 91549.184 516507.090
[53,] 134943.450 91549.184
[54,] 113988.102 134943.450
[55,] 118476.826 113988.102
[56,] -52844.317 118476.826
[57,] 79241.286 -52844.317
[58,] 54817.906 79241.286
[59,] -36863.199 54817.906
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -136858.130 -408463.718
2 -15181.566 -136858.130
3 -295596.702 -15181.566
4 -121363.082 -295596.702
5 3449.126 -121363.082
6 -64776.050 3449.126
7 33905.086 -64776.050
8 83870.631 33905.086
9 69527.901 83870.631
10 6825.466 69527.901
11 -61498.490 6825.466
12 -68551.553 -61498.490
13 -22348.013 -68551.553
14 -157442.345 -22348.013
15 -65120.179 -157442.345
16 41744.328 -65120.179
17 144442.219 41744.328
18 87057.579 144442.219
19 -40990.047 87057.579
20 21564.890 -40990.047
21 -226704.177 21564.890
22 15662.846 -226704.177
23 -76286.701 15662.846
24 59520.212 -76286.701
25 -26760.725 59520.212
26 -32348.488 -26760.725
27 138141.472 -32348.488
28 -94811.896 138141.472
29 198365.702 -94811.896
30 -22710.995 198365.702
31 -49722.035 -22710.995
32 -115287.849 -49722.035
33 -60828.694 -115287.849
34 -100387.377 -60828.694
35 -178991.931 -100387.377
36 24698.039 -178991.931
37 -236377.416 24698.039
38 84092.884 -236377.416
39 285067.493 84092.884
40 -39711.547 285067.493
41 -73813.051 -39711.547
42 252819.023 -73813.051
43 9569.712 252819.023
44 -65022.190 9569.712
45 -31653.038 -65022.190
46 -91248.216 -31653.038
47 -353773.621 -91248.216
48 96254.228 -353773.621
49 316370.944 96254.228
50 341863.712 316370.944
51 516507.090 341863.712
52 91549.184 516507.090
53 134943.450 91549.184
54 113988.102 134943.450
55 118476.826 113988.102
56 -52844.317 118476.826
57 79241.286 -52844.317
58 54817.906 79241.286
59 -36863.199 54817.906
> 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/wessaorg/rcomp/tmp/7qoca1324475533.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/wessaorg/rcomp/tmp/8pqyc1324475533.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/wessaorg/rcomp/tmp/9mezu1324475533.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/wessaorg/rcomp/tmp/10476e1324475533.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11hi4s1324475533.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/wessaorg/rcomp/tmp/12zyap1324475533.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/wessaorg/rcomp/tmp/13ue9d1324475533.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/wessaorg/rcomp/tmp/14bqc21324475533.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/wessaorg/rcomp/tmp/15eoc31324475533.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/wessaorg/rcomp/tmp/16fbf21324475533.tab")
+ }
>
> try(system("convert tmp/1nwfn1324475533.ps tmp/1nwfn1324475533.png",intern=TRUE))
character(0)
> try(system("convert tmp/28q8c1324475533.ps tmp/28q8c1324475533.png",intern=TRUE))
character(0)
> try(system("convert tmp/399cq1324475533.ps tmp/399cq1324475533.png",intern=TRUE))
character(0)
> try(system("convert tmp/4agti1324475533.ps tmp/4agti1324475533.png",intern=TRUE))
character(0)
> try(system("convert tmp/5pv3k1324475533.ps tmp/5pv3k1324475533.png",intern=TRUE))
character(0)
> try(system("convert tmp/6jxgz1324475533.ps tmp/6jxgz1324475533.png",intern=TRUE))
character(0)
> try(system("convert tmp/7qoca1324475533.ps tmp/7qoca1324475533.png",intern=TRUE))
character(0)
> try(system("convert tmp/8pqyc1324475533.ps tmp/8pqyc1324475533.png",intern=TRUE))
character(0)
> try(system("convert tmp/9mezu1324475533.ps tmp/9mezu1324475533.png",intern=TRUE))
character(0)
> try(system("convert tmp/10476e1324475533.ps tmp/10476e1324475533.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.370 0.549 3.932