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(68900
+ ,5960
+ ,44967
+ ,1873
+ ,48500
+ ,9000
+ ,27860
+ ,928
+ ,55500
+ ,9500
+ ,31439
+ ,1126
+ ,62000
+ ,10000
+ ,39592
+ ,1265
+ ,116500
+ ,18000
+ ,72827
+ ,2214
+ ,45000
+ ,8500
+ ,27317
+ ,912
+ ,38000
+ ,8000
+ ,29856
+ ,899
+ ,83000
+ ,23000
+ ,47752
+ ,1803
+ ,59000
+ ,8100
+ ,39117
+ ,1204
+ ,47500
+ ,9000
+ ,29349
+ ,1725
+ ,40500
+ ,7300
+ ,40166
+ ,1080
+ ,40000
+ ,8000
+ ,31679
+ ,1529
+ ,97000
+ ,20000
+ ,58510
+ ,2455
+ ,45500
+ ,8000
+ ,23454
+ ,1151
+ ,40900
+ ,8000
+ ,20897
+ ,1173
+ ,80000
+ ,10500
+ ,56248
+ ,1960
+ ,56000
+ ,4000
+ ,20859
+ ,1344
+ ,37000
+ ,45000
+ ,22610
+ ,988
+ ,50000
+ ,3400
+ ,35947
+ ,1076
+ ,22400
+ ,1500
+ ,5779
+ ,962
+ ,241100
+ ,17800
+ ,50300
+ ,2100
+ ,82200
+ ,18500
+ ,36700
+ ,2300
+ ,234400
+ ,6700
+ ,49100
+ ,1600
+ ,233700
+ ,44200
+ ,52100
+ ,1300
+ ,177700
+ ,3400
+ ,65900
+ ,1700
+ ,65900
+ ,29400
+ ,13800
+ ,2300
+ ,117600
+ ,43200
+ ,28700
+ ,2400
+ ,22500
+ ,2900
+ ,45700
+ ,1000
+ ,326600
+ ,28900
+ ,28400
+ ,1800
+ ,377900
+ ,21000
+ ,7100
+ ,1700
+ ,290700
+ ,8700
+ ,34200
+ ,1400
+ ,108200
+ ,34100
+ ,44200
+ ,1200
+ ,488100
+ ,28100
+ ,5900
+ ,2100
+ ,496600
+ ,38400
+ ,68400
+ ,2200
+ ,493100
+ ,33900
+ ,43600
+ ,1100
+ ,236900
+ ,5100
+ ,66600
+ ,1900
+ ,420600
+ ,14600
+ ,9100
+ ,1500
+ ,328200
+ ,17400
+ ,25500
+ ,2300
+ ,313200
+ ,30900
+ ,8400
+ ,2100
+ ,40200
+ ,25500
+ ,36200
+ ,1800
+ ,318300
+ ,33900
+ ,45000
+ ,1100
+ ,374100
+ ,33700
+ ,11100
+ ,2000
+ ,144400
+ ,19900
+ ,12300
+ ,900
+ ,298300
+ ,26800
+ ,52800
+ ,1500
+ ,404200
+ ,30500
+ ,26000
+ ,1600
+ ,134600
+ ,19500
+ ,9000
+ ,1200
+ ,270600
+ ,42500
+ ,46700
+ ,1000
+ ,181800
+ ,12200
+ ,58200
+ ,2000
+ ,492300
+ ,6600
+ ,54100
+ ,2000
+ ,203000
+ ,2800
+ ,25500
+ ,1900
+ ,464300
+ ,7600
+ ,64500
+ ,2000
+ ,137200
+ ,37700
+ ,42100
+ ,1500
+ ,95100
+ ,28200
+ ,23500
+ ,1900
+ ,481300
+ ,20600
+ ,14000
+ ,1300
+ ,112300
+ ,23000
+ ,65900
+ ,1400
+ ,29500
+ ,15900
+ ,19200
+ ,1300
+ ,76200
+ ,20800
+ ,51200
+ ,1800
+ ,323800
+ ,10000
+ ,50400
+ ,2200
+ ,40600
+ ,22000
+ ,52100
+ ,1000
+ ,425700
+ ,40300
+ ,43400
+ ,1900)
+ ,dim=c(4
+ ,60)
+ ,dimnames=list(c('Verkoopprijs'
+ ,'Grond'
+ ,'waarde'
+ ,'oppervalke')
+ ,1:60))
> y <- array(NA,dim=c(4,60),dimnames=list(c('Verkoopprijs','Grond','waarde','oppervalke'),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
Verkoopprijs Grond waarde oppervalke
1 68900 5960 44967 1873
2 48500 9000 27860 928
3 55500 9500 31439 1126
4 62000 10000 39592 1265
5 116500 18000 72827 2214
6 45000 8500 27317 912
7 38000 8000 29856 899
8 83000 23000 47752 1803
9 59000 8100 39117 1204
10 47500 9000 29349 1725
11 40500 7300 40166 1080
12 40000 8000 31679 1529
13 97000 20000 58510 2455
14 45500 8000 23454 1151
15 40900 8000 20897 1173
16 80000 10500 56248 1960
17 56000 4000 20859 1344
18 37000 45000 22610 988
19 50000 3400 35947 1076
20 22400 1500 5779 962
21 241100 17800 50300 2100
22 82200 18500 36700 2300
23 234400 6700 49100 1600
24 233700 44200 52100 1300
25 177700 3400 65900 1700
26 65900 29400 13800 2300
27 117600 43200 28700 2400
28 22500 2900 45700 1000
29 326600 28900 28400 1800
30 377900 21000 7100 1700
31 290700 8700 34200 1400
32 108200 34100 44200 1200
33 488100 28100 5900 2100
34 496600 38400 68400 2200
35 493100 33900 43600 1100
36 236900 5100 66600 1900
37 420600 14600 9100 1500
38 328200 17400 25500 2300
39 313200 30900 8400 2100
40 40200 25500 36200 1800
41 318300 33900 45000 1100
42 374100 33700 11100 2000
43 144400 19900 12300 900
44 298300 26800 52800 1500
45 404200 30500 26000 1600
46 134600 19500 9000 1200
47 270600 42500 46700 1000
48 181800 12200 58200 2000
49 492300 6600 54100 2000
50 203000 2800 25500 1900
51 464300 7600 64500 2000
52 137200 37700 42100 1500
53 95100 28200 23500 1900
54 481300 20600 14000 1300
55 112300 23000 65900 1400
56 29500 15900 19200 1300
57 76200 20800 51200 1800
58 323800 10000 50400 2200
59 40600 22000 52100 1000
60 425700 40300 43400 1900
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Grond waarde oppervalke
-2.849e+04 3.556e+00 -2.345e-01 9.934e+01
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-239235 -100602 -40805 87902 311322
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.849e+04 7.358e+04 -0.387 0.7001
Grond 3.556e+00 1.507e+00 2.359 0.0218 *
waarde -2.345e-01 1.064e+00 -0.220 0.8264
oppervalke 9.934e+01 4.189e+01 2.371 0.0212 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 142700 on 56 degrees of freedom
Multiple R-squared: 0.1936, Adjusted R-squared: 0.1504
F-statistic: 4.48 on 3 and 56 DF, p-value: 0.00688
> 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,] 3.504228e-06 7.008455e-06 0.9999965
[2,] 3.367019e-06 6.734038e-06 0.9999966
[3,] 1.029475e-07 2.058951e-07 0.9999999
[4,] 3.547585e-09 7.095171e-09 1.0000000
[5,] 6.867618e-09 1.373524e-08 1.0000000
[6,] 7.884992e-10 1.576998e-09 1.0000000
[7,] 4.760144e-11 9.520289e-11 1.0000000
[8,] 4.106941e-12 8.213882e-12 1.0000000
[9,] 2.487806e-13 4.975612e-13 1.0000000
[10,] 1.376674e-14 2.753348e-14 1.0000000
[11,] 1.542162e-14 3.084324e-14 1.0000000
[12,] 2.768748e-15 5.537496e-15 1.0000000
[13,] 1.716540e-16 3.433081e-16 1.0000000
[14,] 1.138015e-17 2.276031e-17 1.0000000
[15,] 2.757451e-08 5.514902e-08 1.0000000
[16,] 1.152665e-08 2.305330e-08 1.0000000
[17,] 7.182153e-07 1.436431e-06 0.9999993
[18,] 1.600504e-06 3.201009e-06 0.9999984
[19,] 6.651184e-07 1.330237e-06 0.9999993
[20,] 5.473647e-07 1.094729e-06 0.9999995
[21,] 6.122468e-07 1.224494e-06 0.9999994
[22,] 4.121513e-07 8.243026e-07 0.9999996
[23,] 4.948935e-05 9.897871e-05 0.9999505
[24,] 2.431920e-03 4.863840e-03 0.9975681
[25,] 5.652802e-03 1.130560e-02 0.9943472
[26,] 3.882193e-03 7.764387e-03 0.9961178
[27,] 2.188518e-02 4.377036e-02 0.9781148
[28,] 5.759965e-02 1.151993e-01 0.9424003
[29,] 2.047685e-01 4.095370e-01 0.7952315
[30,] 1.691578e-01 3.383156e-01 0.8308422
[31,] 2.746665e-01 5.493331e-01 0.7253335
[32,] 2.267760e-01 4.535520e-01 0.7732240
[33,] 1.711366e-01 3.422731e-01 0.8288634
[34,] 2.408742e-01 4.817484e-01 0.7591258
[35,] 2.283389e-01 4.566778e-01 0.7716611
[36,] 1.783859e-01 3.567718e-01 0.8216141
[37,] 1.271261e-01 2.542521e-01 0.8728739
[38,] 9.769013e-02 1.953803e-01 0.9023099
[39,] 1.038393e-01 2.076785e-01 0.8961607
[40,] 6.807044e-02 1.361409e-01 0.9319296
[41,] 6.073586e-02 1.214717e-01 0.9392641
[42,] 4.637231e-02 9.274463e-02 0.9536277
[43,] 8.685694e-02 1.737139e-01 0.9131431
[44,] 5.825541e-02 1.165108e-01 0.9417446
[45,] 1.076687e-01 2.153374e-01 0.8923313
[46,] 6.920895e-02 1.384179e-01 0.9307910
[47,] 1.589243e-01 3.178487e-01 0.8410757
> postscript(file="/var/wessaorg/rcomp/tmp/1x6ba1321808233.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/2uc2q1321808233.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/3t4cx1321808233.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/4bf6j1321808233.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/5mulc1321808233.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
-99327.606 -40675.108 -54283.202 -61457.674 -121888.232 -40934.810
7 8 9 10 11 12
-44269.798 -138221.952 -51752.194 -120499.313 -54843.011 -104426.027
13 14 15 16 17 18
-175799.156 -63304.564 -70689.646 -110370.803 -58359.939 -187397.252
19 20 21 22 23 24
-32065.026 -48657.606 9465.118 -174981.430 91629.459 -11930.534
25 26 27 28 29 30
40671.296 -235416.342 -239234.621 -47949.951 80155.199 164489.868
31 32 33 34 35 36
157190.428 -93429.456 209422.278 186013.674 301975.050 74121.705
37 38 39 40 41 42
250287.725 72304.206 25150.590 -192323.909 127503.351 86659.727
43 44 45 46 47 48
15592.717 94847.253 171369.994 -23360.347 59550.825 -18132.519
49 50 51 52 53 54
311321.907 38763.449 282204.311 -107526.747 -159938.296 310666.186
55 56 57 58 59 60
-64632.512 -123199.295 -136091.279 109994.620 -96276.521 132295.748
> postscript(file="/var/wessaorg/rcomp/tmp/62kai1321808233.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 -99327.606 NA
1 -40675.108 -99327.606
2 -54283.202 -40675.108
3 -61457.674 -54283.202
4 -121888.232 -61457.674
5 -40934.810 -121888.232
6 -44269.798 -40934.810
7 -138221.952 -44269.798
8 -51752.194 -138221.952
9 -120499.313 -51752.194
10 -54843.011 -120499.313
11 -104426.027 -54843.011
12 -175799.156 -104426.027
13 -63304.564 -175799.156
14 -70689.646 -63304.564
15 -110370.803 -70689.646
16 -58359.939 -110370.803
17 -187397.252 -58359.939
18 -32065.026 -187397.252
19 -48657.606 -32065.026
20 9465.118 -48657.606
21 -174981.430 9465.118
22 91629.459 -174981.430
23 -11930.534 91629.459
24 40671.296 -11930.534
25 -235416.342 40671.296
26 -239234.621 -235416.342
27 -47949.951 -239234.621
28 80155.199 -47949.951
29 164489.868 80155.199
30 157190.428 164489.868
31 -93429.456 157190.428
32 209422.278 -93429.456
33 186013.674 209422.278
34 301975.050 186013.674
35 74121.705 301975.050
36 250287.725 74121.705
37 72304.206 250287.725
38 25150.590 72304.206
39 -192323.909 25150.590
40 127503.351 -192323.909
41 86659.727 127503.351
42 15592.717 86659.727
43 94847.253 15592.717
44 171369.994 94847.253
45 -23360.347 171369.994
46 59550.825 -23360.347
47 -18132.519 59550.825
48 311321.907 -18132.519
49 38763.449 311321.907
50 282204.311 38763.449
51 -107526.747 282204.311
52 -159938.296 -107526.747
53 310666.186 -159938.296
54 -64632.512 310666.186
55 -123199.295 -64632.512
56 -136091.279 -123199.295
57 109994.620 -136091.279
58 -96276.521 109994.620
59 132295.748 -96276.521
60 NA 132295.748
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -40675.108 -99327.606
[2,] -54283.202 -40675.108
[3,] -61457.674 -54283.202
[4,] -121888.232 -61457.674
[5,] -40934.810 -121888.232
[6,] -44269.798 -40934.810
[7,] -138221.952 -44269.798
[8,] -51752.194 -138221.952
[9,] -120499.313 -51752.194
[10,] -54843.011 -120499.313
[11,] -104426.027 -54843.011
[12,] -175799.156 -104426.027
[13,] -63304.564 -175799.156
[14,] -70689.646 -63304.564
[15,] -110370.803 -70689.646
[16,] -58359.939 -110370.803
[17,] -187397.252 -58359.939
[18,] -32065.026 -187397.252
[19,] -48657.606 -32065.026
[20,] 9465.118 -48657.606
[21,] -174981.430 9465.118
[22,] 91629.459 -174981.430
[23,] -11930.534 91629.459
[24,] 40671.296 -11930.534
[25,] -235416.342 40671.296
[26,] -239234.621 -235416.342
[27,] -47949.951 -239234.621
[28,] 80155.199 -47949.951
[29,] 164489.868 80155.199
[30,] 157190.428 164489.868
[31,] -93429.456 157190.428
[32,] 209422.278 -93429.456
[33,] 186013.674 209422.278
[34,] 301975.050 186013.674
[35,] 74121.705 301975.050
[36,] 250287.725 74121.705
[37,] 72304.206 250287.725
[38,] 25150.590 72304.206
[39,] -192323.909 25150.590
[40,] 127503.351 -192323.909
[41,] 86659.727 127503.351
[42,] 15592.717 86659.727
[43,] 94847.253 15592.717
[44,] 171369.994 94847.253
[45,] -23360.347 171369.994
[46,] 59550.825 -23360.347
[47,] -18132.519 59550.825
[48,] 311321.907 -18132.519
[49,] 38763.449 311321.907
[50,] 282204.311 38763.449
[51,] -107526.747 282204.311
[52,] -159938.296 -107526.747
[53,] 310666.186 -159938.296
[54,] -64632.512 310666.186
[55,] -123199.295 -64632.512
[56,] -136091.279 -123199.295
[57,] 109994.620 -136091.279
[58,] -96276.521 109994.620
[59,] 132295.748 -96276.521
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -40675.108 -99327.606
2 -54283.202 -40675.108
3 -61457.674 -54283.202
4 -121888.232 -61457.674
5 -40934.810 -121888.232
6 -44269.798 -40934.810
7 -138221.952 -44269.798
8 -51752.194 -138221.952
9 -120499.313 -51752.194
10 -54843.011 -120499.313
11 -104426.027 -54843.011
12 -175799.156 -104426.027
13 -63304.564 -175799.156
14 -70689.646 -63304.564
15 -110370.803 -70689.646
16 -58359.939 -110370.803
17 -187397.252 -58359.939
18 -32065.026 -187397.252
19 -48657.606 -32065.026
20 9465.118 -48657.606
21 -174981.430 9465.118
22 91629.459 -174981.430
23 -11930.534 91629.459
24 40671.296 -11930.534
25 -235416.342 40671.296
26 -239234.621 -235416.342
27 -47949.951 -239234.621
28 80155.199 -47949.951
29 164489.868 80155.199
30 157190.428 164489.868
31 -93429.456 157190.428
32 209422.278 -93429.456
33 186013.674 209422.278
34 301975.050 186013.674
35 74121.705 301975.050
36 250287.725 74121.705
37 72304.206 250287.725
38 25150.590 72304.206
39 -192323.909 25150.590
40 127503.351 -192323.909
41 86659.727 127503.351
42 15592.717 86659.727
43 94847.253 15592.717
44 171369.994 94847.253
45 -23360.347 171369.994
46 59550.825 -23360.347
47 -18132.519 59550.825
48 311321.907 -18132.519
49 38763.449 311321.907
50 282204.311 38763.449
51 -107526.747 282204.311
52 -159938.296 -107526.747
53 310666.186 -159938.296
54 -64632.512 310666.186
55 -123199.295 -64632.512
56 -136091.279 -123199.295
57 109994.620 -136091.279
58 -96276.521 109994.620
59 132295.748 -96276.521
> 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/7at2q1321808233.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/812yy1321808233.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/99a841321808233.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/10nbug1321808233.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/11f6dr1321808233.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/12mkb81321808233.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/13ewkg1321808233.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/14gvzr1321808233.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/15p1hq1321808234.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/163gpd1321808234.tab")
+ }
>
> try(system("convert tmp/1x6ba1321808233.ps tmp/1x6ba1321808233.png",intern=TRUE))
character(0)
> try(system("convert tmp/2uc2q1321808233.ps tmp/2uc2q1321808233.png",intern=TRUE))
character(0)
> try(system("convert tmp/3t4cx1321808233.ps tmp/3t4cx1321808233.png",intern=TRUE))
character(0)
> try(system("convert tmp/4bf6j1321808233.ps tmp/4bf6j1321808233.png",intern=TRUE))
character(0)
> try(system("convert tmp/5mulc1321808233.ps tmp/5mulc1321808233.png",intern=TRUE))
character(0)
> try(system("convert tmp/62kai1321808233.ps tmp/62kai1321808233.png",intern=TRUE))
character(0)
> try(system("convert tmp/7at2q1321808233.ps tmp/7at2q1321808233.png",intern=TRUE))
character(0)
> try(system("convert tmp/812yy1321808233.ps tmp/812yy1321808233.png",intern=TRUE))
character(0)
> try(system("convert tmp/99a841321808233.ps tmp/99a841321808233.png",intern=TRUE))
character(0)
> try(system("convert tmp/10nbug1321808233.ps tmp/10nbug1321808233.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.303 0.541 4.000