R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-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
+ ,1979
+ ,530000
+ ,2700
+ ,345
+ ,1977
+ ,389500
+ ,3571
+ ,150
+ ,1969
+ ,305000
+ ,854
+ ,260
+ ,2011
+ ,620000
+ ,1560
+ ,458
+ ,1981
+ ,750000
+ ,3017
+ ,400
+ ,1972
+ ,389000
+ ,436
+ ,201
+ ,2011
+ ,387000
+ ,1098
+ ,233
+ ,1966
+ ,312000
+ ,625
+ ,160
+ ,1953
+ ,375000
+ ,700
+ ,140
+ ,1985
+ ,385000
+ ,639
+ ,155
+ ,1990
+ ,395000
+ ,1246
+ ,300
+ ,1973
+ ,398000
+ ,600
+ ,220
+ ,1997
+ ,449000
+ ,1000
+ ,272
+ ,1982
+ ,451245
+ ,1047
+ ,280
+ ,2011
+ ,511862
+ ,1414
+ ,219
+ ,2011
+ ,324000
+ ,674
+ ,160
+ ,1973
+ ,772000
+ ,6500
+ ,190
+ ,1971
+ ,617000
+ ,3321
+ ,192
+ ,1980
+ ,595000
+ ,2000
+ ,320
+ ,1997
+ ,475000
+ ,1945
+ ,241
+ ,1974
+ ,985000
+ ,7620
+ ,500
+ ,1977
+ ,439000
+ ,1001
+ ,190
+ ,1991
+ ,479000
+ ,1699
+ ,261
+ ,1987
+ ,657160
+ ,1961
+ ,206
+ ,2012
+ ,299000
+ ,1248
+ ,127
+ ,1970
+ ,419000
+ ,1294
+ ,225
+ ,1989
+ ,449000
+ ,1267
+ ,185
+ ,1970
+ ,327000
+ ,998
+ ,215
+ ,1990
+ ,1695000
+ ,5462
+ ,730
+ ,1998
+ ,489000
+ ,1883
+ ,223
+ ,1987
+ ,449000
+ ,1000
+ ,256
+ ,1994
+ ,470000
+ ,663
+ ,281
+ ,2011
+ ,537000
+ ,2240
+ ,298
+ ,1976
+ ,685000
+ ,2580
+ ,362
+ ,2001
+ ,399000
+ ,2755
+ ,250
+ ,1980
+ ,299500
+ ,773
+ ,188
+ ,1968
+ ,598000
+ ,1465
+ ,500
+ ,1982
+ ,547000
+ ,2025
+ ,270
+ ,1965
+ ,750000
+ ,2160
+ ,300
+ ,1961
+ ,320000
+ ,983
+ ,130
+ ,1989
+ ,373000
+ ,351
+ ,200
+ ,2011
+ ,825000
+ ,712
+ ,270
+ ,2003
+ ,389000
+ ,1120
+ ,224
+ ,1975
+ ,474000
+ ,2619
+ ,290
+ ,1967
+ ,325000
+ ,1193
+ ,214
+ ,1964
+ ,795000
+ ,1500
+ ,450
+ ,2011
+ ,590000
+ ,8560
+ ,330
+ ,1950
+ ,608000
+ ,2236
+ ,190
+ ,1993
+ ,1300000
+ ,3390
+ ,462
+ ,1993
+ ,1325000
+ ,2935
+ ,473
+ ,2004
+ ,1680000
+ ,3700
+ ,528
+ ,2008
+ ,895000
+ ,3290
+ ,470
+ ,1973
+ ,235000
+ ,1115
+ ,94
+ ,1938
+ ,330000
+ ,1200
+ ,100
+ ,1970
+ ,489000
+ ,2160
+ ,166
+ ,1977
+ ,499000
+ ,2605
+ ,334
+ ,1963
+ ,535000
+ ,2229
+ ,230
+ ,1974
+ ,645000
+ ,2267
+ ,303
+ ,1980
+ ,699000
+ ,5027
+ ,315
+ ,1976)
+ ,dim=c(4
+ ,60)
+ ,dimnames=list(c('Verkoopprijs'
+ ,'Oppervlakte'
+ ,'Bewoonbareoppervlakte'
+ ,'Bouwjaar')
+ ,1:60))
> y <- array(NA,dim=c(4,60),dimnames=list(c('Verkoopprijs','Oppervlakte','Bewoonbareoppervlakte','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 = 'Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> 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
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, 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
Verkoopprijs Oppervlakte Bewoonbareoppervlakte Bouwjaar t
1 465000 1520 510 1979 1
2 530000 2700 345 1977 2
3 389500 3571 150 1969 3
4 305000 854 260 2011 4
5 620000 1560 458 1981 5
6 750000 3017 400 1972 6
7 389000 436 201 2011 7
8 387000 1098 233 1966 8
9 312000 625 160 1953 9
10 375000 700 140 1985 10
11 385000 639 155 1990 11
12 395000 1246 300 1973 12
13 398000 600 220 1997 13
14 449000 1000 272 1982 14
15 451245 1047 280 2011 15
16 511862 1414 219 2011 16
17 324000 674 160 1973 17
18 772000 6500 190 1971 18
19 617000 3321 192 1980 19
20 595000 2000 320 1997 20
21 475000 1945 241 1974 21
22 985000 7620 500 1977 22
23 439000 1001 190 1991 23
24 479000 1699 261 1987 24
25 657160 1961 206 2012 25
26 299000 1248 127 1970 26
27 419000 1294 225 1989 27
28 449000 1267 185 1970 28
29 327000 998 215 1990 29
30 1695000 5462 730 1998 30
31 489000 1883 223 1987 31
32 449000 1000 256 1994 32
33 470000 663 281 2011 33
34 537000 2240 298 1976 34
35 685000 2580 362 2001 35
36 399000 2755 250 1980 36
37 299500 773 188 1968 37
38 598000 1465 500 1982 38
39 547000 2025 270 1965 39
40 750000 2160 300 1961 40
41 320000 983 130 1989 41
42 373000 351 200 2011 42
43 825000 712 270 2003 43
44 389000 1120 224 1975 44
45 474000 2619 290 1967 45
46 325000 1193 214 1964 46
47 795000 1500 450 2011 47
48 590000 8560 330 1950 48
49 608000 2236 190 1993 49
50 1300000 3390 462 1993 50
51 1325000 2935 473 2004 51
52 1680000 3700 528 2008 52
53 895000 3290 470 1973 53
54 235000 1115 94 1938 54
55 330000 1200 100 1970 55
56 489000 2160 166 1977 56
57 499000 2605 334 1963 57
58 535000 2229 230 1974 58
59 645000 2267 303 1980 59
60 699000 5027 315 1976 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Oppervlakte Bewoonbareoppervlakte
-9.357e+06 5.804e+01 1.351e+03
Bouwjaar t
4.701e+03 3.455e+03
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-328379 -90425 -24943 92915 489646
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -9.357e+06 2.592e+06 -3.611 0.000661 ***
Oppervlakte 5.804e+01 1.486e+01 3.904 0.000260 ***
Bewoonbareoppervlakte 1.351e+03 1.958e+02 6.901 5.47e-09 ***
Bouwjaar 4.701e+03 1.309e+03 3.590 0.000705 ***
t 3.455e+03 1.184e+03 2.918 0.005099 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 155100 on 55 degrees of freedom
Multiple R-squared: 0.7527, Adjusted R-squared: 0.7347
F-statistic: 41.85 on 4 and 55 DF, p-value: 4.462e-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,] 2.111279e-02 4.222558e-02 0.97888721
[2,] 6.516496e-03 1.303299e-02 0.99348350
[3,] 1.690529e-03 3.381057e-03 0.99830947
[4,] 1.151616e-03 2.303232e-03 0.99884838
[5,] 3.656547e-02 7.313095e-02 0.96343453
[6,] 1.962531e-02 3.925061e-02 0.98037469
[7,] 1.028477e-02 2.056955e-02 0.98971523
[8,] 6.418327e-03 1.283665e-02 0.99358167
[9,] 2.752004e-03 5.504007e-03 0.99724800
[10,] 1.618082e-03 3.236165e-03 0.99838192
[11,] 9.361102e-04 1.872220e-03 0.99906389
[12,] 5.205259e-04 1.041052e-03 0.99947947
[13,] 2.090545e-04 4.181090e-04 0.99979095
[14,] 1.218337e-04 2.436673e-04 0.99987817
[15,] 7.570946e-05 1.514189e-04 0.99992429
[16,] 2.936676e-05 5.873351e-05 0.99997063
[17,] 1.404733e-05 2.809466e-05 0.99998595
[18,] 1.726630e-05 3.453259e-05 0.99998273
[19,] 2.521380e-05 5.042761e-05 0.99997479
[20,] 1.371702e-05 2.743405e-05 0.99998628
[21,] 9.361584e-06 1.872317e-05 0.99999064
[22,] 1.335666e-05 2.671331e-05 0.99998664
[23,] 1.357306e-02 2.714612e-02 0.98642694
[24,] 9.855618e-03 1.971124e-02 0.99014438
[25,] 6.925288e-03 1.385058e-02 0.99307471
[26,] 5.912216e-03 1.182443e-02 0.99408778
[27,] 4.024485e-03 8.048971e-03 0.99597551
[28,] 2.791019e-03 5.582039e-03 0.99720898
[29,] 3.917014e-03 7.834028e-03 0.99608299
[30,] 2.216164e-03 4.432327e-03 0.99778384
[31,] 8.683019e-03 1.736604e-02 0.99131698
[32,] 5.097889e-03 1.019578e-02 0.99490211
[33,] 1.118986e-02 2.237973e-02 0.98881014
[34,] 6.516194e-03 1.303239e-02 0.99348381
[35,] 9.468169e-03 1.893634e-02 0.99053183
[36,] 1.753297e-02 3.506594e-02 0.98246703
[37,] 1.086844e-02 2.173689e-02 0.98913156
[38,] 8.058626e-03 1.611725e-02 0.99194137
[39,] 4.897013e-03 9.794025e-03 0.99510299
[40,] 2.028310e-01 4.056620e-01 0.79716898
[41,] 2.470674e-01 4.941348e-01 0.75293259
[42,] 5.197845e-01 9.604310e-01 0.48021548
[43,] 5.062844e-01 9.874312e-01 0.49371561
[44,] 4.955493e-01 9.910986e-01 0.50445070
[45,] 9.163355e-01 1.673290e-01 0.08366448
> postscript(file="/var/wessaorg/rcomp/tmp/1ly5y1355323797.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/2xkx61355323797.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/3ntvr1355323797.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/4jmbn1355323797.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/50ruj1355323797.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
-262006.297 -36576.116 70031.104 -206321.462 -62280.136 100391.329
7 8 9 10 11 12
-28699.821 95723.796 204475.976 136267.020 102578.762 -42126.827
13 14 15 16 17 18
-9807.556 14769.255 -136303.506 -18010.682 91978.351 167274.004
19 20 21 22 23 24
148302.629 -53368.661 41241.312 -145651.688 42118.853 -38983.028
25 26 27 28 29 30
77317.009 61270.555 -46597.469 124883.368 -119516.089 252434.379
31 32 33 34 35 36
-12494.095 -82202.625 -158796.596 -45214.275 -124405.634 -173952.580
37 38 39 40 41 42
-21690.590 -254223.081 49536.348 219511.517 -47539.073 -159325.507
43 44 45 46 47 48
211285.211 -58063.622 -115091.655 -67986.792 -159105.606 -328379.171
49 50 51 52 53 54
40224.930 294244.345 275621.914 489645.540 -32108.913 103281.882
55 56 57 58 59 60
31358.785 9098.668 -171387.726 -28195.472 -50705.452 -157749.040
> postscript(file="/var/wessaorg/rcomp/tmp/6c9w21355323797.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 -262006.297 NA
1 -36576.116 -262006.297
2 70031.104 -36576.116
3 -206321.462 70031.104
4 -62280.136 -206321.462
5 100391.329 -62280.136
6 -28699.821 100391.329
7 95723.796 -28699.821
8 204475.976 95723.796
9 136267.020 204475.976
10 102578.762 136267.020
11 -42126.827 102578.762
12 -9807.556 -42126.827
13 14769.255 -9807.556
14 -136303.506 14769.255
15 -18010.682 -136303.506
16 91978.351 -18010.682
17 167274.004 91978.351
18 148302.629 167274.004
19 -53368.661 148302.629
20 41241.312 -53368.661
21 -145651.688 41241.312
22 42118.853 -145651.688
23 -38983.028 42118.853
24 77317.009 -38983.028
25 61270.555 77317.009
26 -46597.469 61270.555
27 124883.368 -46597.469
28 -119516.089 124883.368
29 252434.379 -119516.089
30 -12494.095 252434.379
31 -82202.625 -12494.095
32 -158796.596 -82202.625
33 -45214.275 -158796.596
34 -124405.634 -45214.275
35 -173952.580 -124405.634
36 -21690.590 -173952.580
37 -254223.081 -21690.590
38 49536.348 -254223.081
39 219511.517 49536.348
40 -47539.073 219511.517
41 -159325.507 -47539.073
42 211285.211 -159325.507
43 -58063.622 211285.211
44 -115091.655 -58063.622
45 -67986.792 -115091.655
46 -159105.606 -67986.792
47 -328379.171 -159105.606
48 40224.930 -328379.171
49 294244.345 40224.930
50 275621.914 294244.345
51 489645.540 275621.914
52 -32108.913 489645.540
53 103281.882 -32108.913
54 31358.785 103281.882
55 9098.668 31358.785
56 -171387.726 9098.668
57 -28195.472 -171387.726
58 -50705.452 -28195.472
59 -157749.040 -50705.452
60 NA -157749.040
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -36576.116 -262006.297
[2,] 70031.104 -36576.116
[3,] -206321.462 70031.104
[4,] -62280.136 -206321.462
[5,] 100391.329 -62280.136
[6,] -28699.821 100391.329
[7,] 95723.796 -28699.821
[8,] 204475.976 95723.796
[9,] 136267.020 204475.976
[10,] 102578.762 136267.020
[11,] -42126.827 102578.762
[12,] -9807.556 -42126.827
[13,] 14769.255 -9807.556
[14,] -136303.506 14769.255
[15,] -18010.682 -136303.506
[16,] 91978.351 -18010.682
[17,] 167274.004 91978.351
[18,] 148302.629 167274.004
[19,] -53368.661 148302.629
[20,] 41241.312 -53368.661
[21,] -145651.688 41241.312
[22,] 42118.853 -145651.688
[23,] -38983.028 42118.853
[24,] 77317.009 -38983.028
[25,] 61270.555 77317.009
[26,] -46597.469 61270.555
[27,] 124883.368 -46597.469
[28,] -119516.089 124883.368
[29,] 252434.379 -119516.089
[30,] -12494.095 252434.379
[31,] -82202.625 -12494.095
[32,] -158796.596 -82202.625
[33,] -45214.275 -158796.596
[34,] -124405.634 -45214.275
[35,] -173952.580 -124405.634
[36,] -21690.590 -173952.580
[37,] -254223.081 -21690.590
[38,] 49536.348 -254223.081
[39,] 219511.517 49536.348
[40,] -47539.073 219511.517
[41,] -159325.507 -47539.073
[42,] 211285.211 -159325.507
[43,] -58063.622 211285.211
[44,] -115091.655 -58063.622
[45,] -67986.792 -115091.655
[46,] -159105.606 -67986.792
[47,] -328379.171 -159105.606
[48,] 40224.930 -328379.171
[49,] 294244.345 40224.930
[50,] 275621.914 294244.345
[51,] 489645.540 275621.914
[52,] -32108.913 489645.540
[53,] 103281.882 -32108.913
[54,] 31358.785 103281.882
[55,] 9098.668 31358.785
[56,] -171387.726 9098.668
[57,] -28195.472 -171387.726
[58,] -50705.452 -28195.472
[59,] -157749.040 -50705.452
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -36576.116 -262006.297
2 70031.104 -36576.116
3 -206321.462 70031.104
4 -62280.136 -206321.462
5 100391.329 -62280.136
6 -28699.821 100391.329
7 95723.796 -28699.821
8 204475.976 95723.796
9 136267.020 204475.976
10 102578.762 136267.020
11 -42126.827 102578.762
12 -9807.556 -42126.827
13 14769.255 -9807.556
14 -136303.506 14769.255
15 -18010.682 -136303.506
16 91978.351 -18010.682
17 167274.004 91978.351
18 148302.629 167274.004
19 -53368.661 148302.629
20 41241.312 -53368.661
21 -145651.688 41241.312
22 42118.853 -145651.688
23 -38983.028 42118.853
24 77317.009 -38983.028
25 61270.555 77317.009
26 -46597.469 61270.555
27 124883.368 -46597.469
28 -119516.089 124883.368
29 252434.379 -119516.089
30 -12494.095 252434.379
31 -82202.625 -12494.095
32 -158796.596 -82202.625
33 -45214.275 -158796.596
34 -124405.634 -45214.275
35 -173952.580 -124405.634
36 -21690.590 -173952.580
37 -254223.081 -21690.590
38 49536.348 -254223.081
39 219511.517 49536.348
40 -47539.073 219511.517
41 -159325.507 -47539.073
42 211285.211 -159325.507
43 -58063.622 211285.211
44 -115091.655 -58063.622
45 -67986.792 -115091.655
46 -159105.606 -67986.792
47 -328379.171 -159105.606
48 40224.930 -328379.171
49 294244.345 40224.930
50 275621.914 294244.345
51 489645.540 275621.914
52 -32108.913 489645.540
53 103281.882 -32108.913
54 31358.785 103281.882
55 9098.668 31358.785
56 -171387.726 9098.668
57 -28195.472 -171387.726
58 -50705.452 -28195.472
59 -157749.040 -50705.452
> 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/7v7561355323797.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/8ud7f1355323797.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/9jxcf1355323797.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/10r2zv1355323797.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/11cim61355323797.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/12yf0h1355323797.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/13f2zb1355323797.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/14kwsm1355323797.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/15w3gn1355323797.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/16y5121355323797.tab")
+ }
>
> try(system("convert tmp/1ly5y1355323797.ps tmp/1ly5y1355323797.png",intern=TRUE))
character(0)
> try(system("convert tmp/2xkx61355323797.ps tmp/2xkx61355323797.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ntvr1355323797.ps tmp/3ntvr1355323797.png",intern=TRUE))
character(0)
> try(system("convert tmp/4jmbn1355323797.ps tmp/4jmbn1355323797.png",intern=TRUE))
character(0)
> try(system("convert tmp/50ruj1355323797.ps tmp/50ruj1355323797.png",intern=TRUE))
character(0)
> try(system("convert tmp/6c9w21355323797.ps tmp/6c9w21355323797.png",intern=TRUE))
character(0)
> try(system("convert tmp/7v7561355323797.ps tmp/7v7561355323797.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ud7f1355323797.ps tmp/8ud7f1355323797.png",intern=TRUE))
character(0)
> try(system("convert tmp/9jxcf1355323797.ps tmp/9jxcf1355323797.png",intern=TRUE))
character(0)
> try(system("convert tmp/10r2zv1355323797.ps tmp/10r2zv1355323797.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
9.590 1.623 11.234