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(102007
+ ,24776
+ ,4
+ ,0
+ ,112007
+ ,19814
+ ,4
+ ,0
+ ,122007
+ ,12738
+ ,4
+ ,0
+ ,012008
+ ,31566
+ ,4
+ ,0
+ ,022008
+ ,30111
+ ,4
+ ,-32
+ ,032008
+ ,30019
+ ,4
+ ,0
+ ,042008
+ ,31934
+ ,4
+ ,0
+ ,052008
+ ,25826
+ ,4
+ ,0
+ ,062008
+ ,26835
+ ,4
+ ,0
+ ,072008
+ ,20205
+ ,4.18
+ ,0
+ ,082008
+ ,17789
+ ,4.25
+ ,0
+ ,092008
+ ,20520
+ ,4.25
+ ,0
+ ,102008
+ ,22518
+ ,3.97
+ ,0
+ ,112008
+ ,15572
+ ,3.42
+ ,0
+ ,122008
+ ,11509
+ ,2.75
+ ,0
+ ,012009
+ ,25447
+ ,2.31
+ ,0
+ ,022009
+ ,24090
+ ,2
+ ,0
+ ,032009
+ ,27786
+ ,1.66
+ ,0
+ ,042009
+ ,26195
+ ,1.31
+ ,0
+ ,052009
+ ,20516
+ ,1.09
+ ,0
+ ,062009
+ ,22759
+ ,1
+ ,0
+ ,072009
+ ,19028
+ ,1
+ ,0
+ ,082009
+ ,16971
+ ,1
+ ,0
+ ,092009
+ ,20036
+ ,1
+ ,0
+ ,102009
+ ,22485
+ ,1
+ ,0
+ ,112009
+ ,18730
+ ,1
+ ,0
+ ,122009
+ ,14538
+ ,1
+ ,0
+ ,012010
+ ,27561
+ ,1
+ ,0
+ ,022010
+ ,25985
+ ,1
+ ,0
+ ,032010
+ ,34670
+ ,1
+ ,0
+ ,042010
+ ,32066
+ ,1
+ ,0
+ ,052010
+ ,27186
+ ,1
+ ,0
+ ,062010
+ ,29586
+ ,1
+ ,0
+ ,072010
+ ,21359
+ ,1
+ ,0
+ ,082010
+ ,21553
+ ,1
+ ,0
+ ,092010
+ ,19573
+ ,1
+ ,20
+ ,102010
+ ,24256
+ ,1
+ ,0
+ ,112010
+ ,22380
+ ,1
+ ,0
+ ,122010
+ ,16167
+ ,1
+ ,0
+ ,012011
+ ,27297
+ ,1
+ ,0
+ ,022011
+ ,28287
+ ,1
+ ,0
+ ,032011
+ ,33474
+ ,1
+ ,0
+ ,042011
+ ,28229
+ ,1.14
+ ,0
+ ,052011
+ ,28785
+ ,1.25
+ ,0
+ ,062011
+ ,25597
+ ,1.25
+ ,0
+ ,072011
+ ,18130
+ ,1.4
+ ,0
+ ,082011
+ ,20198
+ ,1.5
+ ,0
+ ,092011
+ ,22849
+ ,1.5
+ ,0
+ ,102011
+ ,23118
+ ,1.5
+ ,0
+ ,112011
+ ,21925
+ ,1.32
+ ,0
+ ,122011
+ ,20801
+ ,1.11
+ ,0
+ ,012012
+ ,18785
+ ,1
+ ,0
+ ,122012
+ ,20659
+ ,1
+ ,0
+ ,032012
+ ,29367
+ ,1
+ ,0
+ ,042012
+ ,23992
+ ,1
+ ,0
+ ,052012
+ ,20645
+ ,1
+ ,0
+ ,062012
+ ,22356
+ ,1
+ ,0
+ ,072012
+ ,17902
+ ,0.83
+ ,0
+ ,082012
+ ,15879
+ ,0.75
+ ,0
+ ,092012
+ ,16963
+ ,0.75
+ ,0
+ ,102012
+ ,21035
+ ,0.75
+ ,0)
+ ,dim=c(4
+ ,61)
+ ,dimnames=list(c('Data'
+ ,'Inschrijvingen'
+ ,'Rentevoet'
+ ,'verkoopprijsverloop')
+ ,1:61))
> y <- array(NA,dim=c(4,61),dimnames=list(c('Data','Inschrijvingen','Rentevoet','verkoopprijsverloop'),1:61))
> 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 = '2'
> 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
Inschrijvingen Data Rentevoet verkoopprijsverloop
1 24776 102007 4.00 0
2 19814 112007 4.00 0
3 12738 122007 4.00 0
4 31566 12008 4.00 0
5 30111 22008 4.00 -32
6 30019 32008 4.00 0
7 31934 42008 4.00 0
8 25826 52008 4.00 0
9 26835 62008 4.00 0
10 20205 72008 4.18 0
11 17789 82008 4.25 0
12 20520 92008 4.25 0
13 22518 102008 3.97 0
14 15572 112008 3.42 0
15 11509 122008 2.75 0
16 25447 12009 2.31 0
17 24090 22009 2.00 0
18 27786 32009 1.66 0
19 26195 42009 1.31 0
20 20516 52009 1.09 0
21 22759 62009 1.00 0
22 19028 72009 1.00 0
23 16971 82009 1.00 0
24 20036 92009 1.00 0
25 22485 102009 1.00 0
26 18730 112009 1.00 0
27 14538 122009 1.00 0
28 27561 12010 1.00 0
29 25985 22010 1.00 0
30 34670 32010 1.00 0
31 32066 42010 1.00 0
32 27186 52010 1.00 0
33 29586 62010 1.00 0
34 21359 72010 1.00 0
35 21553 82010 1.00 0
36 19573 92010 1.00 20
37 24256 102010 1.00 0
38 22380 112010 1.00 0
39 16167 122010 1.00 0
40 27297 12011 1.00 0
41 28287 22011 1.00 0
42 33474 32011 1.00 0
43 28229 42011 1.14 0
44 28785 52011 1.25 0
45 25597 62011 1.25 0
46 18130 72011 1.40 0
47 20198 82011 1.50 0
48 22849 92011 1.50 0
49 23118 102011 1.50 0
50 21925 112011 1.32 0
51 20801 122011 1.11 0
52 18785 12012 1.00 0
53 20659 122012 1.00 0
54 29367 32012 1.00 0
55 23992 42012 1.00 0
56 20645 52012 1.00 0
57 22356 62012 1.00 0
58 17902 72012 0.83 0
59 15879 82012 0.75 0
60 16963 92012 0.75 0
61 21035 102012 0.75 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Data Rentevoet
29936.6546 -0.1041 197.2175
verkoopprijsverloop
-51.4764
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-10098.4 -3060.1 46.8 2565.5 7868.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.994e+04 1.276e+03 23.468 < 2e-16 ***
Data -1.041e-01 1.432e-02 -7.268 1.14e-09 ***
Rentevoet 1.972e+02 4.023e+02 0.490 0.626
verkoopprijsverloop -5.148e+01 1.058e+02 -0.486 0.629
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3795 on 57 degrees of freedom
Multiple R-squared: 0.5004, Adjusted R-squared: 0.4741
F-statistic: 19.03 on 3 and 57 DF, p-value: 1.132e-08
> 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,] 0.70707840 0.5858432 0.2929216
[2,] 0.58488040 0.8302392 0.4151196
[3,] 0.45794586 0.9158917 0.5420541
[4,] 0.32476622 0.6495324 0.6752338
[5,] 0.22163197 0.4432639 0.7783680
[6,] 0.21510434 0.4302087 0.7848957
[7,] 0.16838083 0.3367617 0.8316192
[8,] 0.32684418 0.6536884 0.6731558
[9,] 0.31998381 0.6399676 0.6800162
[10,] 0.26410477 0.5282095 0.7358952
[11,] 0.22006744 0.4401349 0.7799326
[12,] 0.26853119 0.5370624 0.7314688
[13,] 0.25812497 0.5162499 0.7418750
[14,] 0.21478536 0.4295707 0.7852146
[15,] 0.18430710 0.3686142 0.8156929
[16,] 0.14557031 0.2911406 0.8544297
[17,] 0.12389271 0.2477854 0.8761073
[18,] 0.11765861 0.2353172 0.8823414
[19,] 0.18691527 0.3738305 0.8130847
[20,] 0.16015847 0.3203169 0.8398415
[21,] 0.13169540 0.2633908 0.8683046
[22,] 0.09455456 0.1891091 0.9054454
[23,] 0.06778262 0.1355652 0.9322174
[24,] 0.26069416 0.5213883 0.7393058
[25,] 0.41800911 0.8360182 0.5819909
[26,] 0.38564732 0.7712946 0.6143527
[27,] 0.52259009 0.9548198 0.4774099
[28,] 0.44889138 0.8977828 0.5511086
[29,] 0.37241684 0.7448337 0.6275832
[30,] 0.29959774 0.5991955 0.7004023
[31,] 0.34443775 0.6888755 0.6555622
[32,] 0.34720332 0.6944066 0.6527967
[33,] 0.28561953 0.5712391 0.7143805
[34,] 0.23146378 0.4629276 0.7685362
[35,] 0.18647394 0.3729479 0.8135261
[36,] 0.46267466 0.9253493 0.5373253
[37,] 0.48599377 0.9719875 0.5140062
[38,] 0.59254191 0.8149162 0.4074581
[39,] 0.57244502 0.8551100 0.4275550
[40,] 0.64407654 0.7118469 0.3559235
[41,] 0.63924512 0.7215098 0.3607549
[42,] 0.55891923 0.8821615 0.4410808
[43,] 0.48779429 0.9755886 0.5122057
[44,] 0.42423056 0.8484611 0.5757694
[45,] 0.34832240 0.6966448 0.6516776
[46,] 0.57396425 0.8520715 0.4260358
[47,] 0.43832275 0.8766455 0.5616772
[48,] 0.74426935 0.5114613 0.2557306
> postscript(file="/var/wessaorg/rcomp/tmp/1npmd1353333913.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/2mhtw1353333913.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/3hlnk1353333913.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/4uwe41353333913.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/5gb941353333913.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 = 61
Frequency = 1
1 2 3 4 5 6
4669.27256 748.25966 -5286.75323 2090.49270 29.23603 2625.46691
7 8 9 10 11 12
5581.45402 514.44112 2564.42823 -3060.08382 -4448.90194 -676.91483
13 14 15 16 17 18
2417.29318 -3379.25008 -6269.12723 -3695.10560 -3949.98107 854.06000
19 20 21 22 23 24
373.07323 -4221.55181 -919.81512 -3609.82802 -4625.84091 -519.85380
25 26 27 28 29 30
2970.13330 256.12041 -2894.89248 -1322.64656 -1857.65945 7868.32766
31 32 33 34 35 36
6305.31476 2466.30187 5907.28898 -1278.72392 -43.73681 46.77765
37 38 39 40 41 42
4741.23740 3906.22451 -1265.78838 -1586.54246 444.44465 6672.43175
43 44 45 46 47 48
2440.80841 4016.10159 1869.08869 -4586.50683 -1497.24147 2194.74564
49 50 51 52 53 54
3504.73274 3388.21900 3346.62179 -10098.43836 3226.41981 2565.53585
55 56 57 58 59 60
-1768.47704 -4074.48993 -1322.50283 -4701.98874 -5668.22423 -3543.23713
61
1569.74998
> postscript(file="/var/wessaorg/rcomp/tmp/6utz51353333913.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 4669.27256 NA
1 748.25966 4669.27256
2 -5286.75323 748.25966
3 2090.49270 -5286.75323
4 29.23603 2090.49270
5 2625.46691 29.23603
6 5581.45402 2625.46691
7 514.44112 5581.45402
8 2564.42823 514.44112
9 -3060.08382 2564.42823
10 -4448.90194 -3060.08382
11 -676.91483 -4448.90194
12 2417.29318 -676.91483
13 -3379.25008 2417.29318
14 -6269.12723 -3379.25008
15 -3695.10560 -6269.12723
16 -3949.98107 -3695.10560
17 854.06000 -3949.98107
18 373.07323 854.06000
19 -4221.55181 373.07323
20 -919.81512 -4221.55181
21 -3609.82802 -919.81512
22 -4625.84091 -3609.82802
23 -519.85380 -4625.84091
24 2970.13330 -519.85380
25 256.12041 2970.13330
26 -2894.89248 256.12041
27 -1322.64656 -2894.89248
28 -1857.65945 -1322.64656
29 7868.32766 -1857.65945
30 6305.31476 7868.32766
31 2466.30187 6305.31476
32 5907.28898 2466.30187
33 -1278.72392 5907.28898
34 -43.73681 -1278.72392
35 46.77765 -43.73681
36 4741.23740 46.77765
37 3906.22451 4741.23740
38 -1265.78838 3906.22451
39 -1586.54246 -1265.78838
40 444.44465 -1586.54246
41 6672.43175 444.44465
42 2440.80841 6672.43175
43 4016.10159 2440.80841
44 1869.08869 4016.10159
45 -4586.50683 1869.08869
46 -1497.24147 -4586.50683
47 2194.74564 -1497.24147
48 3504.73274 2194.74564
49 3388.21900 3504.73274
50 3346.62179 3388.21900
51 -10098.43836 3346.62179
52 3226.41981 -10098.43836
53 2565.53585 3226.41981
54 -1768.47704 2565.53585
55 -4074.48993 -1768.47704
56 -1322.50283 -4074.48993
57 -4701.98874 -1322.50283
58 -5668.22423 -4701.98874
59 -3543.23713 -5668.22423
60 1569.74998 -3543.23713
61 NA 1569.74998
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 748.25966 4669.27256
[2,] -5286.75323 748.25966
[3,] 2090.49270 -5286.75323
[4,] 29.23603 2090.49270
[5,] 2625.46691 29.23603
[6,] 5581.45402 2625.46691
[7,] 514.44112 5581.45402
[8,] 2564.42823 514.44112
[9,] -3060.08382 2564.42823
[10,] -4448.90194 -3060.08382
[11,] -676.91483 -4448.90194
[12,] 2417.29318 -676.91483
[13,] -3379.25008 2417.29318
[14,] -6269.12723 -3379.25008
[15,] -3695.10560 -6269.12723
[16,] -3949.98107 -3695.10560
[17,] 854.06000 -3949.98107
[18,] 373.07323 854.06000
[19,] -4221.55181 373.07323
[20,] -919.81512 -4221.55181
[21,] -3609.82802 -919.81512
[22,] -4625.84091 -3609.82802
[23,] -519.85380 -4625.84091
[24,] 2970.13330 -519.85380
[25,] 256.12041 2970.13330
[26,] -2894.89248 256.12041
[27,] -1322.64656 -2894.89248
[28,] -1857.65945 -1322.64656
[29,] 7868.32766 -1857.65945
[30,] 6305.31476 7868.32766
[31,] 2466.30187 6305.31476
[32,] 5907.28898 2466.30187
[33,] -1278.72392 5907.28898
[34,] -43.73681 -1278.72392
[35,] 46.77765 -43.73681
[36,] 4741.23740 46.77765
[37,] 3906.22451 4741.23740
[38,] -1265.78838 3906.22451
[39,] -1586.54246 -1265.78838
[40,] 444.44465 -1586.54246
[41,] 6672.43175 444.44465
[42,] 2440.80841 6672.43175
[43,] 4016.10159 2440.80841
[44,] 1869.08869 4016.10159
[45,] -4586.50683 1869.08869
[46,] -1497.24147 -4586.50683
[47,] 2194.74564 -1497.24147
[48,] 3504.73274 2194.74564
[49,] 3388.21900 3504.73274
[50,] 3346.62179 3388.21900
[51,] -10098.43836 3346.62179
[52,] 3226.41981 -10098.43836
[53,] 2565.53585 3226.41981
[54,] -1768.47704 2565.53585
[55,] -4074.48993 -1768.47704
[56,] -1322.50283 -4074.48993
[57,] -4701.98874 -1322.50283
[58,] -5668.22423 -4701.98874
[59,] -3543.23713 -5668.22423
[60,] 1569.74998 -3543.23713
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 748.25966 4669.27256
2 -5286.75323 748.25966
3 2090.49270 -5286.75323
4 29.23603 2090.49270
5 2625.46691 29.23603
6 5581.45402 2625.46691
7 514.44112 5581.45402
8 2564.42823 514.44112
9 -3060.08382 2564.42823
10 -4448.90194 -3060.08382
11 -676.91483 -4448.90194
12 2417.29318 -676.91483
13 -3379.25008 2417.29318
14 -6269.12723 -3379.25008
15 -3695.10560 -6269.12723
16 -3949.98107 -3695.10560
17 854.06000 -3949.98107
18 373.07323 854.06000
19 -4221.55181 373.07323
20 -919.81512 -4221.55181
21 -3609.82802 -919.81512
22 -4625.84091 -3609.82802
23 -519.85380 -4625.84091
24 2970.13330 -519.85380
25 256.12041 2970.13330
26 -2894.89248 256.12041
27 -1322.64656 -2894.89248
28 -1857.65945 -1322.64656
29 7868.32766 -1857.65945
30 6305.31476 7868.32766
31 2466.30187 6305.31476
32 5907.28898 2466.30187
33 -1278.72392 5907.28898
34 -43.73681 -1278.72392
35 46.77765 -43.73681
36 4741.23740 46.77765
37 3906.22451 4741.23740
38 -1265.78838 3906.22451
39 -1586.54246 -1265.78838
40 444.44465 -1586.54246
41 6672.43175 444.44465
42 2440.80841 6672.43175
43 4016.10159 2440.80841
44 1869.08869 4016.10159
45 -4586.50683 1869.08869
46 -1497.24147 -4586.50683
47 2194.74564 -1497.24147
48 3504.73274 2194.74564
49 3388.21900 3504.73274
50 3346.62179 3388.21900
51 -10098.43836 3346.62179
52 3226.41981 -10098.43836
53 2565.53585 3226.41981
54 -1768.47704 2565.53585
55 -4074.48993 -1768.47704
56 -1322.50283 -4074.48993
57 -4701.98874 -1322.50283
58 -5668.22423 -4701.98874
59 -3543.23713 -5668.22423
60 1569.74998 -3543.23713
> 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/7fy7i1353333913.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/8ofct1353333913.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/9c3px1353333913.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/10qboo1353333913.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/111cym1353333913.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/1288zm1353333913.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/137n8b1353333913.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/14hgx41353333913.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/15kavv1353333913.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/16mh361353333913.tab")
+ }
>
> try(system("convert tmp/1npmd1353333913.ps tmp/1npmd1353333913.png",intern=TRUE))
character(0)
> try(system("convert tmp/2mhtw1353333913.ps tmp/2mhtw1353333913.png",intern=TRUE))
character(0)
> try(system("convert tmp/3hlnk1353333913.ps tmp/3hlnk1353333913.png",intern=TRUE))
character(0)
> try(system("convert tmp/4uwe41353333913.ps tmp/4uwe41353333913.png",intern=TRUE))
character(0)
> try(system("convert tmp/5gb941353333913.ps tmp/5gb941353333913.png",intern=TRUE))
character(0)
> try(system("convert tmp/6utz51353333913.ps tmp/6utz51353333913.png",intern=TRUE))
character(0)
> try(system("convert tmp/7fy7i1353333913.ps tmp/7fy7i1353333913.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ofct1353333913.ps tmp/8ofct1353333913.png",intern=TRUE))
character(0)
> try(system("convert tmp/9c3px1353333913.ps tmp/9c3px1353333913.png",intern=TRUE))
character(0)
> try(system("convert tmp/10qboo1353333913.ps tmp/10qboo1353333913.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.514 1.058 7.628