R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
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(2173,2363,2126,1905,2121,1983,1734,2074,2049,2406,2558,2251,2059,2397,1747,1707,2319,1631,1627,1791,2034,1997,2169,2028,2253,2218,1855,2187,1852,1570,1851,1954,1828,2251,2277,2085,2282,2266,1878,2267,2069,1746,2299,2360,2214,2825,2355,2333,3016,2155,2172,2150,2533,2058,2160,2259,2498,2695,2799,2945,2930,2318,2540,2570,2669,2450,2842,3439,2677,2979,2257,2842,2546,2455,2293,2379,2478,2054,2272,2351,2271,2542,2304,2194,2722,2395,2146,1894,2548,2087,2063,2481,2476,2212,2834,2148,2598),dim=c(1,97),dimnames=list(c('Bouwvergunningen'),1:97))
> y <- array(NA,dim=c(1,97),dimnames=list(c('Bouwvergunningen'),1:97))
> 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 = 'Include Monthly 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.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
Bouwvergunningen M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 2173 1 0 0 0 0 0 0 0 0 0 0
2 2363 0 1 0 0 0 0 0 0 0 0 0
3 2126 0 0 1 0 0 0 0 0 0 0 0
4 1905 0 0 0 1 0 0 0 0 0 0 0
5 2121 0 0 0 0 1 0 0 0 0 0 0
6 1983 0 0 0 0 0 1 0 0 0 0 0
7 1734 0 0 0 0 0 0 1 0 0 0 0
8 2074 0 0 0 0 0 0 0 1 0 0 0
9 2049 0 0 0 0 0 0 0 0 1 0 0
10 2406 0 0 0 0 0 0 0 0 0 1 0
11 2558 0 0 0 0 0 0 0 0 0 0 1
12 2251 0 0 0 0 0 0 0 0 0 0 0
13 2059 1 0 0 0 0 0 0 0 0 0 0
14 2397 0 1 0 0 0 0 0 0 0 0 0
15 1747 0 0 1 0 0 0 0 0 0 0 0
16 1707 0 0 0 1 0 0 0 0 0 0 0
17 2319 0 0 0 0 1 0 0 0 0 0 0
18 1631 0 0 0 0 0 1 0 0 0 0 0
19 1627 0 0 0 0 0 0 1 0 0 0 0
20 1791 0 0 0 0 0 0 0 1 0 0 0
21 2034 0 0 0 0 0 0 0 0 1 0 0
22 1997 0 0 0 0 0 0 0 0 0 1 0
23 2169 0 0 0 0 0 0 0 0 0 0 1
24 2028 0 0 0 0 0 0 0 0 0 0 0
25 2253 1 0 0 0 0 0 0 0 0 0 0
26 2218 0 1 0 0 0 0 0 0 0 0 0
27 1855 0 0 1 0 0 0 0 0 0 0 0
28 2187 0 0 0 1 0 0 0 0 0 0 0
29 1852 0 0 0 0 1 0 0 0 0 0 0
30 1570 0 0 0 0 0 1 0 0 0 0 0
31 1851 0 0 0 0 0 0 1 0 0 0 0
32 1954 0 0 0 0 0 0 0 1 0 0 0
33 1828 0 0 0 0 0 0 0 0 1 0 0
34 2251 0 0 0 0 0 0 0 0 0 1 0
35 2277 0 0 0 0 0 0 0 0 0 0 1
36 2085 0 0 0 0 0 0 0 0 0 0 0
37 2282 1 0 0 0 0 0 0 0 0 0 0
38 2266 0 1 0 0 0 0 0 0 0 0 0
39 1878 0 0 1 0 0 0 0 0 0 0 0
40 2267 0 0 0 1 0 0 0 0 0 0 0
41 2069 0 0 0 0 1 0 0 0 0 0 0
42 1746 0 0 0 0 0 1 0 0 0 0 0
43 2299 0 0 0 0 0 0 1 0 0 0 0
44 2360 0 0 0 0 0 0 0 1 0 0 0
45 2214 0 0 0 0 0 0 0 0 1 0 0
46 2825 0 0 0 0 0 0 0 0 0 1 0
47 2355 0 0 0 0 0 0 0 0 0 0 1
48 2333 0 0 0 0 0 0 0 0 0 0 0
49 3016 1 0 0 0 0 0 0 0 0 0 0
50 2155 0 1 0 0 0 0 0 0 0 0 0
51 2172 0 0 1 0 0 0 0 0 0 0 0
52 2150 0 0 0 1 0 0 0 0 0 0 0
53 2533 0 0 0 0 1 0 0 0 0 0 0
54 2058 0 0 0 0 0 1 0 0 0 0 0
55 2160 0 0 0 0 0 0 1 0 0 0 0
56 2259 0 0 0 0 0 0 0 1 0 0 0
57 2498 0 0 0 0 0 0 0 0 1 0 0
58 2695 0 0 0 0 0 0 0 0 0 1 0
59 2799 0 0 0 0 0 0 0 0 0 0 1
60 2945 0 0 0 0 0 0 0 0 0 0 0
61 2930 1 0 0 0 0 0 0 0 0 0 0
62 2318 0 1 0 0 0 0 0 0 0 0 0
63 2540 0 0 1 0 0 0 0 0 0 0 0
64 2570 0 0 0 1 0 0 0 0 0 0 0
65 2669 0 0 0 0 1 0 0 0 0 0 0
66 2450 0 0 0 0 0 1 0 0 0 0 0
67 2842 0 0 0 0 0 0 1 0 0 0 0
68 3439 0 0 0 0 0 0 0 1 0 0 0
69 2677 0 0 0 0 0 0 0 0 1 0 0
70 2979 0 0 0 0 0 0 0 0 0 1 0
71 2257 0 0 0 0 0 0 0 0 0 0 1
72 2842 0 0 0 0 0 0 0 0 0 0 0
73 2546 1 0 0 0 0 0 0 0 0 0 0
74 2455 0 1 0 0 0 0 0 0 0 0 0
75 2293 0 0 1 0 0 0 0 0 0 0 0
76 2379 0 0 0 1 0 0 0 0 0 0 0
77 2478 0 0 0 0 1 0 0 0 0 0 0
78 2054 0 0 0 0 0 1 0 0 0 0 0
79 2272 0 0 0 0 0 0 1 0 0 0 0
80 2351 0 0 0 0 0 0 0 1 0 0 0
81 2271 0 0 0 0 0 0 0 0 1 0 0
82 2542 0 0 0 0 0 0 0 0 0 1 0
83 2304 0 0 0 0 0 0 0 0 0 0 1
84 2194 0 0 0 0 0 0 0 0 0 0 0
85 2722 1 0 0 0 0 0 0 0 0 0 0
86 2395 0 1 0 0 0 0 0 0 0 0 0
87 2146 0 0 1 0 0 0 0 0 0 0 0
88 1894 0 0 0 1 0 0 0 0 0 0 0
89 2548 0 0 0 0 1 0 0 0 0 0 0
90 2087 0 0 0 0 0 1 0 0 0 0 0
91 2063 0 0 0 0 0 0 1 0 0 0 0
92 2481 0 0 0 0 0 0 0 1 0 0 0
93 2476 0 0 0 0 0 0 0 0 1 0 0
94 2212 0 0 0 0 0 0 0 0 0 1 0
95 2834 0 0 0 0 0 0 0 0 0 0 1
96 2148 0 0 0 0 0 0 0 0 0 0 0
97 2598 1 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
2353.25 155.53 -32.38 -258.63 -220.87 -29.63
M6 M7 M8 M9 M10 M11
-405.87 -247.25 -14.63 -97.38 135.12 90.87
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-547.63 -227.37 12.37 166.00 1100.38
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2353.25 112.59 20.902 <2e-16 ***
M1 155.53 154.73 1.005 0.3177
M2 -32.38 159.22 -0.203 0.8394
M3 -258.63 159.22 -1.624 0.1080
M4 -220.87 159.22 -1.387 0.1690
M5 -29.63 159.22 -0.186 0.8528
M6 -405.87 159.22 -2.549 0.0126 *
M7 -247.25 159.22 -1.553 0.1242
M8 -14.63 159.22 -0.092 0.9270
M9 -97.38 159.22 -0.612 0.5425
M10 135.12 159.22 0.849 0.3985
M11 90.87 159.22 0.571 0.5697
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 318.4 on 85 degrees of freedom
Multiple R-squared: 0.2395, Adjusted R-squared: 0.141
F-statistic: 2.433 on 11 and 85 DF, p-value: 0.01094
> 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.132547601 0.265095202 0.867452399
[2,] 0.079230197 0.158460394 0.920769803
[3,] 0.046007857 0.092015715 0.953992143
[4,] 0.053746580 0.107493160 0.946253420
[5,] 0.030219701 0.060439401 0.969780299
[6,] 0.029735567 0.059471135 0.970264433
[7,] 0.014500251 0.029000502 0.985499749
[8,] 0.025689592 0.051379183 0.974310408
[9,] 0.032024395 0.064048789 0.967975605
[10,] 0.024410243 0.048820487 0.975589757
[11,] 0.016850444 0.033700887 0.983149556
[12,] 0.010904075 0.021808149 0.989095925
[13,] 0.006537142 0.013074284 0.993462858
[14,] 0.009802852 0.019605704 0.990197148
[15,] 0.016588156 0.033176311 0.983411844
[16,] 0.016330247 0.032660494 0.983669753
[17,] 0.013791643 0.027583286 0.986208357
[18,] 0.011973132 0.023946265 0.988026868
[19,] 0.013597565 0.027195129 0.986402435
[20,] 0.009933154 0.019866308 0.990066846
[21,] 0.006524046 0.013048091 0.993475954
[22,] 0.004824666 0.009649332 0.995175334
[23,] 0.004350324 0.008700647 0.995649676
[24,] 0.002551220 0.005102441 0.997448780
[25,] 0.001891215 0.003782430 0.998108785
[26,] 0.002486616 0.004973232 0.997513384
[27,] 0.002168906 0.004337811 0.997831094
[28,] 0.001713060 0.003426120 0.998286940
[29,] 0.007578462 0.015156923 0.992421538
[30,] 0.012882989 0.025765978 0.987117011
[31,] 0.012168762 0.024337524 0.987831238
[32,] 0.034683393 0.069366786 0.965316607
[33,] 0.025082633 0.050165266 0.974917367
[34,] 0.020895062 0.041790124 0.979104938
[35,] 0.106786424 0.213572848 0.893213576
[36,] 0.089463252 0.178926504 0.910536748
[37,] 0.078862078 0.157724156 0.921137922
[38,] 0.060784051 0.121568103 0.939215949
[39,] 0.065757887 0.131515773 0.934242113
[40,] 0.059931491 0.119862982 0.940068509
[41,] 0.054205229 0.108410458 0.945794771
[42,] 0.062273623 0.124547245 0.937726377
[43,] 0.066080060 0.132160119 0.933919940
[44,] 0.057899045 0.115798090 0.942100955
[45,] 0.069792309 0.139584618 0.930207691
[46,] 0.171306199 0.342612399 0.828693801
[47,] 0.207621994 0.415243987 0.792378006
[48,] 0.163264249 0.326528498 0.836735751
[49,] 0.193915025 0.387830049 0.806084975
[50,] 0.233454255 0.466908511 0.766545745
[51,] 0.223247581 0.446495162 0.776752419
[52,] 0.270330799 0.540661598 0.729669201
[53,] 0.488837906 0.977675812 0.511162094
[54,] 0.941195125 0.117609750 0.058804875
[55,] 0.940219841 0.119560319 0.059780159
[56,] 0.975585302 0.048829396 0.024414698
[57,] 0.972355504 0.055288992 0.027644496
[58,] 0.996763125 0.006473750 0.003236875
[59,] 0.993741717 0.012516566 0.006258283
[60,] 0.987404553 0.025190894 0.012595447
[61,] 0.977886073 0.044227854 0.022113927
[62,] 0.988367193 0.023265613 0.011632807
[63,] 0.975668784 0.048662432 0.024331216
[64,] 0.949973525 0.100052951 0.050026475
[65,] 0.920421149 0.159157702 0.079578851
[66,] 0.858096481 0.283807038 0.141903519
[67,] 0.778272767 0.443454466 0.221727233
[68,] 0.726943867 0.546112265 0.273056133
> postscript(file="/var/www/html/rcomp/tmp/1bfzx1227473645.ps",horizontal=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/html/rcomp/tmp/26zud1227473645.ps",horizontal=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/html/rcomp/tmp/3527g1227473645.ps",horizontal=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/html/rcomp/tmp/40glc1227473645.ps",horizontal=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/html/rcomp/tmp/5mqvr1227473645.ps",horizontal=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 = 97
Frequency = 1
1 2 3 4 5 6 7
-335.77778 42.12500 31.37500 -227.37500 -202.62500 35.62500 -372.00000
8 9 10 11 12 13 14
-264.62500 -206.87500 -82.37500 113.87500 -102.25000 -449.77778 76.12500
15 16 17 18 19 20 21
-347.62500 -425.37500 -4.62500 -316.37500 -479.00000 -547.62500 -221.87500
22 23 24 25 26 27 28
-491.37500 -275.12500 -325.25000 -255.77778 -102.87500 -239.62500 54.62500
29 30 31 32 33 34 35
-471.62500 -377.37500 -255.00000 -384.62500 -427.87500 -237.37500 -167.12500
36 37 38 39 40 41 42
-268.25000 -226.77778 -54.87500 -216.62500 134.62500 -254.62500 -201.37500
43 44 45 46 47 48 49
193.00000 21.37500 -41.87500 336.62500 -89.12500 -20.25000 507.22222
50 51 52 53 54 55 56
-165.87500 77.37500 17.62500 209.37500 110.62500 54.00000 -79.62500
57 58 59 60 61 62 63
242.12500 206.62500 354.87500 591.75000 421.22222 -2.87500 445.37500
64 65 66 67 68 69 70
437.62500 345.37500 502.62500 736.00000 1100.37500 421.12500 490.62500
71 72 73 74 75 76 77
-187.12500 488.75000 37.22222 134.12500 198.37500 246.62500 154.37500
78 79 80 81 82 83 84
106.62500 166.00000 12.37500 15.12500 53.62500 -140.12500 -159.25000
85 86 87 88 89 90 91
213.22222 74.12500 51.37500 -238.37500 224.37500 139.62500 -43.00000
92 93 94 95 96 97
142.37500 220.12500 -276.37500 389.87500 -205.25000 89.22222
> postscript(file="/var/www/html/rcomp/tmp/6tw801227473645.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 97
Frequency = 1
lag(myerror, k = 1) myerror
0 -335.77778 NA
1 42.12500 -335.77778
2 31.37500 42.12500
3 -227.37500 31.37500
4 -202.62500 -227.37500
5 35.62500 -202.62500
6 -372.00000 35.62500
7 -264.62500 -372.00000
8 -206.87500 -264.62500
9 -82.37500 -206.87500
10 113.87500 -82.37500
11 -102.25000 113.87500
12 -449.77778 -102.25000
13 76.12500 -449.77778
14 -347.62500 76.12500
15 -425.37500 -347.62500
16 -4.62500 -425.37500
17 -316.37500 -4.62500
18 -479.00000 -316.37500
19 -547.62500 -479.00000
20 -221.87500 -547.62500
21 -491.37500 -221.87500
22 -275.12500 -491.37500
23 -325.25000 -275.12500
24 -255.77778 -325.25000
25 -102.87500 -255.77778
26 -239.62500 -102.87500
27 54.62500 -239.62500
28 -471.62500 54.62500
29 -377.37500 -471.62500
30 -255.00000 -377.37500
31 -384.62500 -255.00000
32 -427.87500 -384.62500
33 -237.37500 -427.87500
34 -167.12500 -237.37500
35 -268.25000 -167.12500
36 -226.77778 -268.25000
37 -54.87500 -226.77778
38 -216.62500 -54.87500
39 134.62500 -216.62500
40 -254.62500 134.62500
41 -201.37500 -254.62500
42 193.00000 -201.37500
43 21.37500 193.00000
44 -41.87500 21.37500
45 336.62500 -41.87500
46 -89.12500 336.62500
47 -20.25000 -89.12500
48 507.22222 -20.25000
49 -165.87500 507.22222
50 77.37500 -165.87500
51 17.62500 77.37500
52 209.37500 17.62500
53 110.62500 209.37500
54 54.00000 110.62500
55 -79.62500 54.00000
56 242.12500 -79.62500
57 206.62500 242.12500
58 354.87500 206.62500
59 591.75000 354.87500
60 421.22222 591.75000
61 -2.87500 421.22222
62 445.37500 -2.87500
63 437.62500 445.37500
64 345.37500 437.62500
65 502.62500 345.37500
66 736.00000 502.62500
67 1100.37500 736.00000
68 421.12500 1100.37500
69 490.62500 421.12500
70 -187.12500 490.62500
71 488.75000 -187.12500
72 37.22222 488.75000
73 134.12500 37.22222
74 198.37500 134.12500
75 246.62500 198.37500
76 154.37500 246.62500
77 106.62500 154.37500
78 166.00000 106.62500
79 12.37500 166.00000
80 15.12500 12.37500
81 53.62500 15.12500
82 -140.12500 53.62500
83 -159.25000 -140.12500
84 213.22222 -159.25000
85 74.12500 213.22222
86 51.37500 74.12500
87 -238.37500 51.37500
88 224.37500 -238.37500
89 139.62500 224.37500
90 -43.00000 139.62500
91 142.37500 -43.00000
92 220.12500 142.37500
93 -276.37500 220.12500
94 389.87500 -276.37500
95 -205.25000 389.87500
96 89.22222 -205.25000
97 NA 89.22222
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 42.12500 -335.77778
[2,] 31.37500 42.12500
[3,] -227.37500 31.37500
[4,] -202.62500 -227.37500
[5,] 35.62500 -202.62500
[6,] -372.00000 35.62500
[7,] -264.62500 -372.00000
[8,] -206.87500 -264.62500
[9,] -82.37500 -206.87500
[10,] 113.87500 -82.37500
[11,] -102.25000 113.87500
[12,] -449.77778 -102.25000
[13,] 76.12500 -449.77778
[14,] -347.62500 76.12500
[15,] -425.37500 -347.62500
[16,] -4.62500 -425.37500
[17,] -316.37500 -4.62500
[18,] -479.00000 -316.37500
[19,] -547.62500 -479.00000
[20,] -221.87500 -547.62500
[21,] -491.37500 -221.87500
[22,] -275.12500 -491.37500
[23,] -325.25000 -275.12500
[24,] -255.77778 -325.25000
[25,] -102.87500 -255.77778
[26,] -239.62500 -102.87500
[27,] 54.62500 -239.62500
[28,] -471.62500 54.62500
[29,] -377.37500 -471.62500
[30,] -255.00000 -377.37500
[31,] -384.62500 -255.00000
[32,] -427.87500 -384.62500
[33,] -237.37500 -427.87500
[34,] -167.12500 -237.37500
[35,] -268.25000 -167.12500
[36,] -226.77778 -268.25000
[37,] -54.87500 -226.77778
[38,] -216.62500 -54.87500
[39,] 134.62500 -216.62500
[40,] -254.62500 134.62500
[41,] -201.37500 -254.62500
[42,] 193.00000 -201.37500
[43,] 21.37500 193.00000
[44,] -41.87500 21.37500
[45,] 336.62500 -41.87500
[46,] -89.12500 336.62500
[47,] -20.25000 -89.12500
[48,] 507.22222 -20.25000
[49,] -165.87500 507.22222
[50,] 77.37500 -165.87500
[51,] 17.62500 77.37500
[52,] 209.37500 17.62500
[53,] 110.62500 209.37500
[54,] 54.00000 110.62500
[55,] -79.62500 54.00000
[56,] 242.12500 -79.62500
[57,] 206.62500 242.12500
[58,] 354.87500 206.62500
[59,] 591.75000 354.87500
[60,] 421.22222 591.75000
[61,] -2.87500 421.22222
[62,] 445.37500 -2.87500
[63,] 437.62500 445.37500
[64,] 345.37500 437.62500
[65,] 502.62500 345.37500
[66,] 736.00000 502.62500
[67,] 1100.37500 736.00000
[68,] 421.12500 1100.37500
[69,] 490.62500 421.12500
[70,] -187.12500 490.62500
[71,] 488.75000 -187.12500
[72,] 37.22222 488.75000
[73,] 134.12500 37.22222
[74,] 198.37500 134.12500
[75,] 246.62500 198.37500
[76,] 154.37500 246.62500
[77,] 106.62500 154.37500
[78,] 166.00000 106.62500
[79,] 12.37500 166.00000
[80,] 15.12500 12.37500
[81,] 53.62500 15.12500
[82,] -140.12500 53.62500
[83,] -159.25000 -140.12500
[84,] 213.22222 -159.25000
[85,] 74.12500 213.22222
[86,] 51.37500 74.12500
[87,] -238.37500 51.37500
[88,] 224.37500 -238.37500
[89,] 139.62500 224.37500
[90,] -43.00000 139.62500
[91,] 142.37500 -43.00000
[92,] 220.12500 142.37500
[93,] -276.37500 220.12500
[94,] 389.87500 -276.37500
[95,] -205.25000 389.87500
[96,] 89.22222 -205.25000
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 42.12500 -335.77778
2 31.37500 42.12500
3 -227.37500 31.37500
4 -202.62500 -227.37500
5 35.62500 -202.62500
6 -372.00000 35.62500
7 -264.62500 -372.00000
8 -206.87500 -264.62500
9 -82.37500 -206.87500
10 113.87500 -82.37500
11 -102.25000 113.87500
12 -449.77778 -102.25000
13 76.12500 -449.77778
14 -347.62500 76.12500
15 -425.37500 -347.62500
16 -4.62500 -425.37500
17 -316.37500 -4.62500
18 -479.00000 -316.37500
19 -547.62500 -479.00000
20 -221.87500 -547.62500
21 -491.37500 -221.87500
22 -275.12500 -491.37500
23 -325.25000 -275.12500
24 -255.77778 -325.25000
25 -102.87500 -255.77778
26 -239.62500 -102.87500
27 54.62500 -239.62500
28 -471.62500 54.62500
29 -377.37500 -471.62500
30 -255.00000 -377.37500
31 -384.62500 -255.00000
32 -427.87500 -384.62500
33 -237.37500 -427.87500
34 -167.12500 -237.37500
35 -268.25000 -167.12500
36 -226.77778 -268.25000
37 -54.87500 -226.77778
38 -216.62500 -54.87500
39 134.62500 -216.62500
40 -254.62500 134.62500
41 -201.37500 -254.62500
42 193.00000 -201.37500
43 21.37500 193.00000
44 -41.87500 21.37500
45 336.62500 -41.87500
46 -89.12500 336.62500
47 -20.25000 -89.12500
48 507.22222 -20.25000
49 -165.87500 507.22222
50 77.37500 -165.87500
51 17.62500 77.37500
52 209.37500 17.62500
53 110.62500 209.37500
54 54.00000 110.62500
55 -79.62500 54.00000
56 242.12500 -79.62500
57 206.62500 242.12500
58 354.87500 206.62500
59 591.75000 354.87500
60 421.22222 591.75000
61 -2.87500 421.22222
62 445.37500 -2.87500
63 437.62500 445.37500
64 345.37500 437.62500
65 502.62500 345.37500
66 736.00000 502.62500
67 1100.37500 736.00000
68 421.12500 1100.37500
69 490.62500 421.12500
70 -187.12500 490.62500
71 488.75000 -187.12500
72 37.22222 488.75000
73 134.12500 37.22222
74 198.37500 134.12500
75 246.62500 198.37500
76 154.37500 246.62500
77 106.62500 154.37500
78 166.00000 106.62500
79 12.37500 166.00000
80 15.12500 12.37500
81 53.62500 15.12500
82 -140.12500 53.62500
83 -159.25000 -140.12500
84 213.22222 -159.25000
85 74.12500 213.22222
86 51.37500 74.12500
87 -238.37500 51.37500
88 224.37500 -238.37500
89 139.62500 224.37500
90 -43.00000 139.62500
91 142.37500 -43.00000
92 220.12500 142.37500
93 -276.37500 220.12500
94 389.87500 -276.37500
95 -205.25000 389.87500
96 89.22222 -205.25000
> 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/html/rcomp/tmp/7wxn21227473645.ps",horizontal=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/html/rcomp/tmp/8sox51227473645.ps",horizontal=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/html/rcomp/tmp/9dklc1227473645.ps",horizontal=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/html/rcomp/tmp/10x8j61227473645.ps",horizontal=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/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/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/html/rcomp/tmp/116kgf1227473645.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/html/rcomp/tmp/12b2ro1227473645.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/html/rcomp/tmp/13ferk1227473645.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/html/rcomp/tmp/14ddnh1227473645.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/html/rcomp/tmp/15lb261227473645.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/html/rcomp/tmp/161z3r1227473645.tab")
+ }
>
> system("convert tmp/1bfzx1227473645.ps tmp/1bfzx1227473645.png")
> system("convert tmp/26zud1227473645.ps tmp/26zud1227473645.png")
> system("convert tmp/3527g1227473645.ps tmp/3527g1227473645.png")
> system("convert tmp/40glc1227473645.ps tmp/40glc1227473645.png")
> system("convert tmp/5mqvr1227473645.ps tmp/5mqvr1227473645.png")
> system("convert tmp/6tw801227473645.ps tmp/6tw801227473645.png")
> system("convert tmp/7wxn21227473645.ps tmp/7wxn21227473645.png")
> system("convert tmp/8sox51227473645.ps tmp/8sox51227473645.png")
> system("convert tmp/9dklc1227473645.ps tmp/9dklc1227473645.png")
> system("convert tmp/10x8j61227473645.ps tmp/10x8j61227473645.png")
>
>
> proc.time()
user system elapsed
2.818 1.594 3.289