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(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'
+ ,'Grondwaarde'
+ ,'Waardehuis'
+ ,'Oppervlakte')
+ ,1:60))
> y <- array(NA,dim=c(4,60),dimnames=list(c('Verkoopprijs','Grondwaarde','Waardehuis','Oppervlakte'),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'
> #'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
Verkoopprijs Grondwaarde Waardehuis Oppervlakte t
1 68900 5960 44967 1873 1
2 48500 9000 27860 928 2
3 55500 9500 31439 1126 3
4 62000 10000 39592 1265 4
5 116500 18000 72827 2214 5
6 45000 8500 27317 912 6
7 38000 8000 29856 899 7
8 83000 23000 47752 1803 8
9 59000 8100 39117 1204 9
10 47500 9000 29349 1725 10
11 40500 7300 40166 1080 11
12 40000 8000 31679 1529 12
13 97000 20000 58510 2455 13
14 45500 8000 23454 1151 14
15 40900 8000 20897 1173 15
16 80000 10500 56248 1960 16
17 56000 4000 20859 1344 17
18 37000 45000 22610 988 18
19 50000 3400 35947 1076 19
20 22400 1500 5779 962 20
21 241100 17800 50300 2100 21
22 82200 18500 36700 2300 22
23 234400 6700 49100 1600 23
24 233700 44200 52100 1300 24
25 177700 3400 65900 1700 25
26 65900 29400 13800 2300 26
27 117600 43200 28700 2400 27
28 22500 2900 45700 1000 28
29 326600 28900 28400 1800 29
30 377900 21000 7100 1700 30
31 290700 8700 34200 1400 31
32 108200 34100 44200 1200 32
33 488100 28100 5900 2100 33
34 496600 38400 68400 2200 34
35 493100 33900 43600 1100 35
36 236900 5100 66600 1900 36
37 420600 14600 9100 1500 37
38 328200 17400 25500 2300 38
39 313200 30900 8400 2100 39
40 40200 25500 36200 1800 40
41 318300 33900 45000 1100 41
42 374100 33700 11100 2000 42
43 144400 19900 12300 900 43
44 298300 26800 52800 1500 44
45 404200 30500 26000 1600 45
46 134600 19500 9000 1200 46
47 270600 42500 46700 1000 47
48 181800 12200 58200 2000 48
49 492300 6600 54100 2000 49
50 203000 2800 25500 1900 50
51 464300 7600 64500 2000 51
52 137200 37700 42100 1500 52
53 95100 28200 23500 1900 53
54 481300 20600 14000 1300 54
55 112300 23000 65900 1400 55
56 29500 15900 19200 1300 56
57 76200 20800 51200 1800 57
58 323800 10000 50400 2200 58
59 40600 22000 52100 1000 59
60 425700 40300 43400 1900 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Grondwaarde Waardehuis Oppervlakte t
-7.827e+04 1.856e+00 -3.074e-01 8.427e+01 3.554e+03
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-224384 -75265 -5319 62526 304788
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.827e+04 6.920e+04 -1.131 0.26290
Grondwaarde 1.856e+00 1.474e+00 1.259 0.21331
Waardehuis -3.074e-01 9.780e-01 -0.314 0.75446
Oppervlakte 8.427e+01 3.875e+01 2.175 0.03396 *
t 3.554e+03 1.055e+03 3.368 0.00139 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 131100 on 55 degrees of freedom
Multiple R-squared: 0.3315, Adjusted R-squared: 0.2828
F-statistic: 6.817 on 4 and 55 DF, p-value: 0.0001570
> 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.369921e-05 4.739842e-05 0.9999763
[2,] 1.690323e-06 3.380646e-06 0.9999983
[3,] 5.044893e-08 1.008979e-07 0.9999999
[4,] 1.229571e-08 2.459143e-08 1.0000000
[5,] 4.217796e-10 8.435591e-10 1.0000000
[6,] 1.996188e-11 3.992375e-11 1.0000000
[7,] 1.956234e-11 3.912467e-11 1.0000000
[8,] 1.962800e-12 3.925599e-12 1.0000000
[9,] 1.057185e-13 2.114370e-13 1.0000000
[10,] 2.264349e-13 4.528699e-13 1.0000000
[11,] 3.008223e-14 6.016446e-14 1.0000000
[12,] 1.828394e-15 3.656788e-15 1.0000000
[13,] 1.206195e-16 2.412391e-16 1.0000000
[14,] 3.682556e-08 7.365112e-08 1.0000000
[15,] 1.596402e-08 3.192804e-08 1.0000000
[16,] 9.942636e-08 1.988527e-07 0.9999999
[17,] 6.187246e-08 1.237449e-07 0.9999999
[18,] 2.128551e-08 4.257101e-08 1.0000000
[19,] 1.486444e-08 2.972887e-08 1.0000000
[20,] 1.996855e-08 3.993709e-08 1.0000000
[21,] 2.857811e-07 5.715622e-07 0.9999997
[22,] 2.005444e-05 4.010887e-05 0.9999799
[23,] 7.301916e-04 1.460383e-03 0.9992698
[24,] 6.551331e-04 1.310266e-03 0.9993449
[25,] 1.006900e-03 2.013799e-03 0.9989931
[26,] 5.798354e-03 1.159671e-02 0.9942016
[27,] 1.029122e-02 2.058244e-02 0.9897088
[28,] 2.809618e-02 5.619236e-02 0.9719038
[29,] 2.028416e-02 4.056832e-02 0.9797158
[30,] 2.408910e-02 4.817819e-02 0.9759109
[31,] 1.469634e-02 2.939267e-02 0.9853037
[32,] 8.848988e-03 1.769798e-02 0.9911510
[33,] 7.435524e-02 1.487105e-01 0.9256448
[34,] 5.055781e-02 1.011156e-01 0.9494422
[35,] 3.232311e-02 6.464622e-02 0.9676769
[36,] 2.524391e-02 5.048783e-02 0.9747561
[37,] 1.470578e-02 2.941156e-02 0.9852942
[38,] 1.040288e-02 2.080575e-02 0.9895971
[39,] 8.466587e-03 1.693317e-02 0.9915334
[40,] 5.896818e-03 1.179364e-02 0.9941032
[41,] 6.034779e-03 1.206956e-02 0.9939652
[42,] 7.465169e-03 1.493034e-02 0.9925348
[43,] 4.780782e-03 9.561565e-03 0.9952192
[44,] 1.468930e-02 2.937860e-02 0.9853107
[45,] 8.959395e-03 1.791879e-02 0.9910406
> postscript(file="/var/www/rcomp/tmp/1sxwm1321980637.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/2qqo71321980637.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/3bouw1321980637.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/4ny0p1321980637.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/5hhzs1321980637.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
-11458.3645 33323.0753 20256.2198 13067.3390 -20589.6887 17718.2761
7 8 9 10 11 12
9968.6742 -47102.5818 820.3238 -62811.0158 -12530.0556 -58328.9768
13 14 15 16 17 18
-96939.9234 -30610.1757 -41403.6713 -65951.0259 -40409.0822 -108515.8416
19 20 21 22 23 24
-25180.1944 -52474.1900 50207.0038 -134580.3034 98766.6786 51120.3489
25 26 27 28 29 30
37821.2474 -192363.0595 -173674.8161 -74331.5338 105227.4169 169514.9852
31 32 33 34 35 36
135200.4668 -78064.8009 221800.4858 218416.3466 304788.3603 38138.3929
37 38 39 40 41 42
216686.8046 53161.7315 21151.2280 -211553.7939 109097.4762 75451.1046
43 44 45 46 47 48
-39124.7930 60303.2459 139117.6804 -85138.5667 33065.1453 -83790.0909
49 50 51 52 53 54
232289.0673 -53876.5170 198522.9966 -152743.4539 -220191.6731 224201.5685
55 56 57 58 59 60
-145279.3781 -224384.3491 -222630.3358 7505.7330 -199871.5369 69194.3664
> postscript(file="/var/www/rcomp/tmp/6o01b1321980637.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 -11458.3645 NA
1 33323.0753 -11458.3645
2 20256.2198 33323.0753
3 13067.3390 20256.2198
4 -20589.6887 13067.3390
5 17718.2761 -20589.6887
6 9968.6742 17718.2761
7 -47102.5818 9968.6742
8 820.3238 -47102.5818
9 -62811.0158 820.3238
10 -12530.0556 -62811.0158
11 -58328.9768 -12530.0556
12 -96939.9234 -58328.9768
13 -30610.1757 -96939.9234
14 -41403.6713 -30610.1757
15 -65951.0259 -41403.6713
16 -40409.0822 -65951.0259
17 -108515.8416 -40409.0822
18 -25180.1944 -108515.8416
19 -52474.1900 -25180.1944
20 50207.0038 -52474.1900
21 -134580.3034 50207.0038
22 98766.6786 -134580.3034
23 51120.3489 98766.6786
24 37821.2474 51120.3489
25 -192363.0595 37821.2474
26 -173674.8161 -192363.0595
27 -74331.5338 -173674.8161
28 105227.4169 -74331.5338
29 169514.9852 105227.4169
30 135200.4668 169514.9852
31 -78064.8009 135200.4668
32 221800.4858 -78064.8009
33 218416.3466 221800.4858
34 304788.3603 218416.3466
35 38138.3929 304788.3603
36 216686.8046 38138.3929
37 53161.7315 216686.8046
38 21151.2280 53161.7315
39 -211553.7939 21151.2280
40 109097.4762 -211553.7939
41 75451.1046 109097.4762
42 -39124.7930 75451.1046
43 60303.2459 -39124.7930
44 139117.6804 60303.2459
45 -85138.5667 139117.6804
46 33065.1453 -85138.5667
47 -83790.0909 33065.1453
48 232289.0673 -83790.0909
49 -53876.5170 232289.0673
50 198522.9966 -53876.5170
51 -152743.4539 198522.9966
52 -220191.6731 -152743.4539
53 224201.5685 -220191.6731
54 -145279.3781 224201.5685
55 -224384.3491 -145279.3781
56 -222630.3358 -224384.3491
57 7505.7330 -222630.3358
58 -199871.5369 7505.7330
59 69194.3664 -199871.5369
60 NA 69194.3664
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 33323.0753 -11458.3645
[2,] 20256.2198 33323.0753
[3,] 13067.3390 20256.2198
[4,] -20589.6887 13067.3390
[5,] 17718.2761 -20589.6887
[6,] 9968.6742 17718.2761
[7,] -47102.5818 9968.6742
[8,] 820.3238 -47102.5818
[9,] -62811.0158 820.3238
[10,] -12530.0556 -62811.0158
[11,] -58328.9768 -12530.0556
[12,] -96939.9234 -58328.9768
[13,] -30610.1757 -96939.9234
[14,] -41403.6713 -30610.1757
[15,] -65951.0259 -41403.6713
[16,] -40409.0822 -65951.0259
[17,] -108515.8416 -40409.0822
[18,] -25180.1944 -108515.8416
[19,] -52474.1900 -25180.1944
[20,] 50207.0038 -52474.1900
[21,] -134580.3034 50207.0038
[22,] 98766.6786 -134580.3034
[23,] 51120.3489 98766.6786
[24,] 37821.2474 51120.3489
[25,] -192363.0595 37821.2474
[26,] -173674.8161 -192363.0595
[27,] -74331.5338 -173674.8161
[28,] 105227.4169 -74331.5338
[29,] 169514.9852 105227.4169
[30,] 135200.4668 169514.9852
[31,] -78064.8009 135200.4668
[32,] 221800.4858 -78064.8009
[33,] 218416.3466 221800.4858
[34,] 304788.3603 218416.3466
[35,] 38138.3929 304788.3603
[36,] 216686.8046 38138.3929
[37,] 53161.7315 216686.8046
[38,] 21151.2280 53161.7315
[39,] -211553.7939 21151.2280
[40,] 109097.4762 -211553.7939
[41,] 75451.1046 109097.4762
[42,] -39124.7930 75451.1046
[43,] 60303.2459 -39124.7930
[44,] 139117.6804 60303.2459
[45,] -85138.5667 139117.6804
[46,] 33065.1453 -85138.5667
[47,] -83790.0909 33065.1453
[48,] 232289.0673 -83790.0909
[49,] -53876.5170 232289.0673
[50,] 198522.9966 -53876.5170
[51,] -152743.4539 198522.9966
[52,] -220191.6731 -152743.4539
[53,] 224201.5685 -220191.6731
[54,] -145279.3781 224201.5685
[55,] -224384.3491 -145279.3781
[56,] -222630.3358 -224384.3491
[57,] 7505.7330 -222630.3358
[58,] -199871.5369 7505.7330
[59,] 69194.3664 -199871.5369
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 33323.0753 -11458.3645
2 20256.2198 33323.0753
3 13067.3390 20256.2198
4 -20589.6887 13067.3390
5 17718.2761 -20589.6887
6 9968.6742 17718.2761
7 -47102.5818 9968.6742
8 820.3238 -47102.5818
9 -62811.0158 820.3238
10 -12530.0556 -62811.0158
11 -58328.9768 -12530.0556
12 -96939.9234 -58328.9768
13 -30610.1757 -96939.9234
14 -41403.6713 -30610.1757
15 -65951.0259 -41403.6713
16 -40409.0822 -65951.0259
17 -108515.8416 -40409.0822
18 -25180.1944 -108515.8416
19 -52474.1900 -25180.1944
20 50207.0038 -52474.1900
21 -134580.3034 50207.0038
22 98766.6786 -134580.3034
23 51120.3489 98766.6786
24 37821.2474 51120.3489
25 -192363.0595 37821.2474
26 -173674.8161 -192363.0595
27 -74331.5338 -173674.8161
28 105227.4169 -74331.5338
29 169514.9852 105227.4169
30 135200.4668 169514.9852
31 -78064.8009 135200.4668
32 221800.4858 -78064.8009
33 218416.3466 221800.4858
34 304788.3603 218416.3466
35 38138.3929 304788.3603
36 216686.8046 38138.3929
37 53161.7315 216686.8046
38 21151.2280 53161.7315
39 -211553.7939 21151.2280
40 109097.4762 -211553.7939
41 75451.1046 109097.4762
42 -39124.7930 75451.1046
43 60303.2459 -39124.7930
44 139117.6804 60303.2459
45 -85138.5667 139117.6804
46 33065.1453 -85138.5667
47 -83790.0909 33065.1453
48 232289.0673 -83790.0909
49 -53876.5170 232289.0673
50 198522.9966 -53876.5170
51 -152743.4539 198522.9966
52 -220191.6731 -152743.4539
53 224201.5685 -220191.6731
54 -145279.3781 224201.5685
55 -224384.3491 -145279.3781
56 -222630.3358 -224384.3491
57 7505.7330 -222630.3358
58 -199871.5369 7505.7330
59 69194.3664 -199871.5369
> 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/7d2571321980637.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/8i0la1321980637.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/9yhwo1321980637.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/107lvj1321980637.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/119qfw1321980637.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/12ogl51321980637.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/13ylym1321980637.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/14lzyz1321980637.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/15236k1321980637.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/164ru21321980637.tab")
+ }
>
> try(system("convert tmp/1sxwm1321980637.ps tmp/1sxwm1321980637.png",intern=TRUE))
character(0)
> try(system("convert tmp/2qqo71321980637.ps tmp/2qqo71321980637.png",intern=TRUE))
character(0)
> try(system("convert tmp/3bouw1321980637.ps tmp/3bouw1321980637.png",intern=TRUE))
character(0)
> try(system("convert tmp/4ny0p1321980637.ps tmp/4ny0p1321980637.png",intern=TRUE))
character(0)
> try(system("convert tmp/5hhzs1321980637.ps tmp/5hhzs1321980637.png",intern=TRUE))
character(0)
> try(system("convert tmp/6o01b1321980637.ps tmp/6o01b1321980637.png",intern=TRUE))
character(0)
> try(system("convert tmp/7d2571321980637.ps tmp/7d2571321980637.png",intern=TRUE))
character(0)
> try(system("convert tmp/8i0la1321980637.ps tmp/8i0la1321980637.png",intern=TRUE))
character(0)
> try(system("convert tmp/9yhwo1321980637.ps tmp/9yhwo1321980637.png",intern=TRUE))
character(0)
> try(system("convert tmp/107lvj1321980637.ps tmp/107lvj1321980637.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.890 0.430 4.303