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(180144
+ ,966.2
+ ,173666
+ ,1153.2
+ ,165688
+ ,1328.3
+ ,161570
+ ,1144.5
+ ,156145
+ ,1477.1
+ ,153730
+ ,1234.9
+ ,182698
+ ,1119.1
+ ,200765
+ ,1356.9
+ ,176512
+ ,1217
+ ,166618
+ ,1440.5
+ ,158644
+ ,1556.6
+ ,159585
+ ,1303.6
+ ,163095
+ ,1421.5
+ ,159044
+ ,1172.5
+ ,155511
+ ,1422.1
+ ,153745
+ ,1263
+ ,150569
+ ,1428.1
+ ,150605
+ ,1347
+ ,179612
+ ,1224.2
+ ,194690
+ ,1201.3
+ ,189917
+ ,997.8
+ ,184128
+ ,1248.8
+ ,175335
+ ,1268.6
+ ,179566
+ ,1016.7
+ ,181140
+ ,1194.3
+ ,177876
+ ,1181.8
+ ,175041
+ ,1150.7
+ ,169292
+ ,1247.2
+ ,166070
+ ,1260.6
+ ,166972
+ ,1249.3
+ ,206348
+ ,1223.2
+ ,215706
+ ,1153
+ ,202108
+ ,1191.5
+ ,195411
+ ,1303.1
+ ,193111
+ ,1267.1
+ ,195198
+ ,1125.2
+ ,198770
+ ,1322.4
+ ,194163
+ ,1089.2
+ ,190420
+ ,1147.3
+ ,189733
+ ,1196.4
+ ,186029
+ ,1190.2
+ ,191531
+ ,1146
+ ,232571
+ ,1139.8
+ ,243477
+ ,1045.6
+ ,227247
+ ,1050.9
+ ,217859
+ ,1117.3
+ ,208679
+ ,1120
+ ,213188
+ ,1052.1
+ ,216234
+ ,1065.8
+ ,213586
+ ,1092.5
+ ,209465
+ ,1422
+ ,204045
+ ,1367.5
+ ,200237
+ ,1136.3
+ ,203666
+ ,1293.7
+ ,241476
+ ,1154.8
+ ,260307
+ ,1206.7
+ ,243324
+ ,1199
+ ,244460
+ ,1265
+ ,233575
+ ,1247.1
+ ,237217
+ ,1116.5
+ ,235243
+ ,1153.9
+ ,230354
+ ,1077.4
+ ,227184
+ ,1132.5
+ ,221678
+ ,1058.8
+ ,217142
+ ,1195.1
+ ,219452
+ ,1263.4
+ ,256446
+ ,1023.1
+ ,265845
+ ,1141
+ ,248624
+ ,1116.3
+ ,241114
+ ,1135.6
+ ,229245
+ ,1210.5
+ ,231805
+ ,1230
+ ,219277
+ ,1136.5
+ ,219313
+ ,1068.7
+ ,212610
+ ,1372.5
+ ,214771
+ ,1049.9
+ ,211142
+ ,1302.2
+ ,211457
+ ,1305.9
+ ,240048
+ ,1173.5
+ ,240636
+ ,1277.4
+ ,230580
+ ,1238.6
+ ,208795
+ ,1508.6
+ ,197922
+ ,1423.4
+ ,194596
+ ,1375.1)
+ ,dim=c(2
+ ,84)
+ ,dimnames=list(c('werkloosheid'
+ ,'Amerika')
+ ,1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('werkloosheid','Amerika'),1:84))
> 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
werkloosheid Amerika M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 180144 966.2 1 0 0 0 0 0 0 0 0 0 0
2 173666 1153.2 0 1 0 0 0 0 0 0 0 0 0
3 165688 1328.3 0 0 1 0 0 0 0 0 0 0 0
4 161570 1144.5 0 0 0 1 0 0 0 0 0 0 0
5 156145 1477.1 0 0 0 0 1 0 0 0 0 0 0
6 153730 1234.9 0 0 0 0 0 1 0 0 0 0 0
7 182698 1119.1 0 0 0 0 0 0 1 0 0 0 0
8 200765 1356.9 0 0 0 0 0 0 0 1 0 0 0
9 176512 1217.0 0 0 0 0 0 0 0 0 1 0 0
10 166618 1440.5 0 0 0 0 0 0 0 0 0 1 0
11 158644 1556.6 0 0 0 0 0 0 0 0 0 0 1
12 159585 1303.6 0 0 0 0 0 0 0 0 0 0 0
13 163095 1421.5 1 0 0 0 0 0 0 0 0 0 0
14 159044 1172.5 0 1 0 0 0 0 0 0 0 0 0
15 155511 1422.1 0 0 1 0 0 0 0 0 0 0 0
16 153745 1263.0 0 0 0 1 0 0 0 0 0 0 0
17 150569 1428.1 0 0 0 0 1 0 0 0 0 0 0
18 150605 1347.0 0 0 0 0 0 1 0 0 0 0 0
19 179612 1224.2 0 0 0 0 0 0 1 0 0 0 0
20 194690 1201.3 0 0 0 0 0 0 0 1 0 0 0
21 189917 997.8 0 0 0 0 0 0 0 0 1 0 0
22 184128 1248.8 0 0 0 0 0 0 0 0 0 1 0
23 175335 1268.6 0 0 0 0 0 0 0 0 0 0 1
24 179566 1016.7 0 0 0 0 0 0 0 0 0 0 0
25 181140 1194.3 1 0 0 0 0 0 0 0 0 0 0
26 177876 1181.8 0 1 0 0 0 0 0 0 0 0 0
27 175041 1150.7 0 0 1 0 0 0 0 0 0 0 0
28 169292 1247.2 0 0 0 1 0 0 0 0 0 0 0
29 166070 1260.6 0 0 0 0 1 0 0 0 0 0 0
30 166972 1249.3 0 0 0 0 0 1 0 0 0 0 0
31 206348 1223.2 0 0 0 0 0 0 1 0 0 0 0
32 215706 1153.0 0 0 0 0 0 0 0 1 0 0 0
33 202108 1191.5 0 0 0 0 0 0 0 0 1 0 0
34 195411 1303.1 0 0 0 0 0 0 0 0 0 1 0
35 193111 1267.1 0 0 0 0 0 0 0 0 0 0 1
36 195198 1125.2 0 0 0 0 0 0 0 0 0 0 0
37 198770 1322.4 1 0 0 0 0 0 0 0 0 0 0
38 194163 1089.2 0 1 0 0 0 0 0 0 0 0 0
39 190420 1147.3 0 0 1 0 0 0 0 0 0 0 0
40 189733 1196.4 0 0 0 1 0 0 0 0 0 0 0
41 186029 1190.2 0 0 0 0 1 0 0 0 0 0 0
42 191531 1146.0 0 0 0 0 0 1 0 0 0 0 0
43 232571 1139.8 0 0 0 0 0 0 1 0 0 0 0
44 243477 1045.6 0 0 0 0 0 0 0 1 0 0 0
45 227247 1050.9 0 0 0 0 0 0 0 0 1 0 0
46 217859 1117.3 0 0 0 0 0 0 0 0 0 1 0
47 208679 1120.0 0 0 0 0 0 0 0 0 0 0 1
48 213188 1052.1 0 0 0 0 0 0 0 0 0 0 0
49 216234 1065.8 1 0 0 0 0 0 0 0 0 0 0
50 213586 1092.5 0 1 0 0 0 0 0 0 0 0 0
51 209465 1422.0 0 0 1 0 0 0 0 0 0 0 0
52 204045 1367.5 0 0 0 1 0 0 0 0 0 0 0
53 200237 1136.3 0 0 0 0 1 0 0 0 0 0 0
54 203666 1293.7 0 0 0 0 0 1 0 0 0 0 0
55 241476 1154.8 0 0 0 0 0 0 1 0 0 0 0
56 260307 1206.7 0 0 0 0 0 0 0 1 0 0 0
57 243324 1199.0 0 0 0 0 0 0 0 0 1 0 0
58 244460 1265.0 0 0 0 0 0 0 0 0 0 1 0
59 233575 1247.1 0 0 0 0 0 0 0 0 0 0 1
60 237217 1116.5 0 0 0 0 0 0 0 0 0 0 0
61 235243 1153.9 1 0 0 0 0 0 0 0 0 0 0
62 230354 1077.4 0 1 0 0 0 0 0 0 0 0 0
63 227184 1132.5 0 0 1 0 0 0 0 0 0 0 0
64 221678 1058.8 0 0 0 1 0 0 0 0 0 0 0
65 217142 1195.1 0 0 0 0 1 0 0 0 0 0 0
66 219452 1263.4 0 0 0 0 0 1 0 0 0 0 0
67 256446 1023.1 0 0 0 0 0 0 1 0 0 0 0
68 265845 1141.0 0 0 0 0 0 0 0 1 0 0 0
69 248624 1116.3 0 0 0 0 0 0 0 0 1 0 0
70 241114 1135.6 0 0 0 0 0 0 0 0 0 1 0
71 229245 1210.5 0 0 0 0 0 0 0 0 0 0 1
72 231805 1230.0 0 0 0 0 0 0 0 0 0 0 0
73 219277 1136.5 1 0 0 0 0 0 0 0 0 0 0
74 219313 1068.7 0 1 0 0 0 0 0 0 0 0 0
75 212610 1372.5 0 0 1 0 0 0 0 0 0 0 0
76 214771 1049.9 0 0 0 1 0 0 0 0 0 0 0
77 211142 1302.2 0 0 0 0 1 0 0 0 0 0 0
78 211457 1305.9 0 0 0 0 0 1 0 0 0 0 0
79 240048 1173.5 0 0 0 0 0 0 1 0 0 0 0
80 240636 1277.4 0 0 0 0 0 0 0 1 0 0 0
81 230580 1238.6 0 0 0 0 0 0 0 0 1 0 0
82 208795 1508.6 0 0 0 0 0 0 0 0 0 1 0
83 197922 1423.4 0 0 0 0 0 0 0 0 0 0 1
84 194596 1375.1 0 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) Amerika M1 M2 M3 M4
312055.29 -94.08 -1908.18 -11324.13 -585.07 -12307.34
M5 M6 M7 M8 M9 M10
-7334.94 -7902.94 16121.53 32225.32 12511.39 17494.69
M11
9655.44
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-40780 -19865 -3269 22018 35464
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 312055.29 32113.49 9.717 1.13e-14 ***
Amerika -94.08 26.08 -3.607 0.000573 ***
M1 -1908.18 13660.96 -0.140 0.889307
M2 -11324.13 13734.79 -0.824 0.412425
M3 -585.07 13947.68 -0.042 0.966658
M4 -12307.34 13666.03 -0.901 0.370857
M5 -7334.94 13958.47 -0.525 0.600884
M6 -7902.94 13854.69 -0.570 0.570196
M7 16121.53 13673.34 1.179 0.242315
M8 32225.32 13673.54 2.357 0.021197 *
M9 12511.39 13682.08 0.914 0.363583
M10 17494.69 13981.33 1.251 0.214936
M11 9655.44 14043.03 0.688 0.493969
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 25560 on 71 degrees of freedom
Multiple R-squared: 0.3559, Adjusted R-squared: 0.2471
F-statistic: 3.27 on 12 and 71 DF, p-value: 0.0008625
> 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.568675e-02 5.137350e-02 9.743132e-01
[2,] 8.931435e-03 1.786287e-02 9.910686e-01
[3,] 2.277986e-03 4.555973e-03 9.977220e-01
[4,] 6.131456e-04 1.226291e-03 9.993869e-01
[5,] 5.666913e-04 1.133383e-03 9.994333e-01
[6,] 2.295216e-04 4.590432e-04 9.997705e-01
[7,] 1.491428e-04 2.982857e-04 9.998509e-01
[8,] 5.467110e-05 1.093422e-04 9.999453e-01
[9,] 2.828998e-05 5.657996e-05 9.999717e-01
[10,] 2.584026e-05 5.168051e-05 9.999742e-01
[11,] 3.890925e-05 7.781850e-05 9.999611e-01
[12,] 2.153380e-05 4.306759e-05 9.999785e-01
[13,] 4.572406e-05 9.144813e-05 9.999543e-01
[14,] 2.706640e-05 5.413281e-05 9.999729e-01
[15,] 5.909448e-05 1.181890e-04 9.999409e-01
[16,] 1.500459e-03 3.000918e-03 9.984995e-01
[17,] 2.328307e-03 4.656614e-03 9.976717e-01
[18,] 9.810109e-03 1.962022e-02 9.901899e-01
[19,] 1.946828e-02 3.893655e-02 9.805317e-01
[20,] 3.088806e-02 6.177612e-02 9.691119e-01
[21,] 6.813662e-02 1.362732e-01 9.318634e-01
[22,] 1.641191e-01 3.282383e-01 8.358809e-01
[23,] 2.336376e-01 4.672753e-01 7.663624e-01
[24,] 3.365194e-01 6.730388e-01 6.634806e-01
[25,] 4.762296e-01 9.524591e-01 5.237704e-01
[26,] 5.485049e-01 9.029902e-01 4.514951e-01
[27,] 7.014698e-01 5.970604e-01 2.985302e-01
[28,] 8.358330e-01 3.283340e-01 1.641670e-01
[29,] 8.827711e-01 2.344577e-01 1.172289e-01
[30,] 9.334175e-01 1.331649e-01 6.658247e-02
[31,] 9.616321e-01 7.673589e-02 3.836794e-02
[32,] 9.774992e-01 4.500155e-02 2.250078e-02
[33,] 9.917709e-01 1.645827e-02 8.229133e-03
[34,] 9.944817e-01 1.103651e-02 5.518257e-03
[35,] 9.954773e-01 9.045384e-03 4.522692e-03
[36,] 9.980915e-01 3.816930e-03 1.908465e-03
[37,] 9.993633e-01 1.273342e-03 6.366709e-04
[38,] 9.998993e-01 2.014568e-04 1.007284e-04
[39,] 9.999231e-01 1.537758e-04 7.688788e-05
[40,] 9.998797e-01 2.405762e-04 1.202881e-04
[41,] 9.998668e-01 2.663019e-04 1.331510e-04
[42,] 9.998400e-01 3.199829e-04 1.599914e-04
[43,] 9.998774e-01 2.452254e-04 1.226127e-04
[44,] 9.998717e-01 2.566998e-04 1.283499e-04
[45,] 9.997452e-01 5.096971e-04 2.548485e-04
[46,] 9.998080e-01 3.840328e-04 1.920164e-04
[47,] 9.997023e-01 5.953778e-04 2.976889e-04
[48,] 9.995472e-01 9.056155e-04 4.528078e-04
[49,] 9.988228e-01 2.354463e-03 1.177231e-03
[50,] 9.968328e-01 6.334472e-03 3.167236e-03
[51,] 9.905592e-01 1.888170e-02 9.440848e-03
[52,] 9.706714e-01 5.865711e-02 2.932856e-02
[53,] 9.276989e-01 1.446022e-01 7.230112e-02
> postscript(file="/var/www/html/rcomp/tmp/1rcg51229727827.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/2619j1229727827.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/39fg81229727827.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/43quo1229727827.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/57h281229727827.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 = 84
Frequency = 1
1 2 3 4 5 6
-39106.5796 -18576.3580 -20820.6534 -30507.6143 -9615.2241 -34247.5071
7 8 9 10 11 12
-40198.0135 -15863.4576 -33563.8057 -27415.0509 -16627.5448 -29832.4085
13 14 15 16 17 18
-13322.6338 -31382.6851 -22173.2951 -27184.5711 -19800.9634 -26826.5523
19 20 21 22 23 24
-33396.5929 -36576.7320 -40780.3338 -27939.4803 -27030.5232 -36841.9030
25 26 27 28 29 30
-16651.7724 -11675.7754 -28175.6067 -13123.9768 -20057.7460 -19650.8082
31 32 33 34 35 36
-6754.6692 -20104.6180 -10366.7517 -11548.1364 -9395.6377 -11002.6229
37 38 39 40 41 42
13029.4034 -4100.2421 -13116.4662 2537.9464 -6721.7186 -4809.8914
43 44 45 46 47 48
11622.3662 -2437.4141 1545.1185 -6579.5156 -7666.2635 110.3985
49 50 51 52 53 54
6353.4213 15633.2097 31771.2973 32946.4037 2415.5681 21220.1801
55 56 57 58 59 60
21938.5109 29548.2801 31554.8206 33916.5560 29186.8360 30197.9131
61 62 63 64 65 66
33650.5445 30980.6574 22255.2044 21538.0456 24852.2554 34155.6678
67 68 69 70 71 72
24518.6604 28905.4662 29074.7094 18397.0810 21413.6429 35463.5748
73 74 75 76 77 78
16047.6167 19121.1935 30259.5197 13793.7664 28927.8286 30158.9112
79 80 81 82 83 84
22269.7380 16528.4755 22536.2427 21168.5461 10119.4902 11905.0480
> postscript(file="/var/www/html/rcomp/tmp/6xhev1229727827.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 = 84
Frequency = 1
lag(myerror, k = 1) myerror
0 -39106.5796 NA
1 -18576.3580 -39106.5796
2 -20820.6534 -18576.3580
3 -30507.6143 -20820.6534
4 -9615.2241 -30507.6143
5 -34247.5071 -9615.2241
6 -40198.0135 -34247.5071
7 -15863.4576 -40198.0135
8 -33563.8057 -15863.4576
9 -27415.0509 -33563.8057
10 -16627.5448 -27415.0509
11 -29832.4085 -16627.5448
12 -13322.6338 -29832.4085
13 -31382.6851 -13322.6338
14 -22173.2951 -31382.6851
15 -27184.5711 -22173.2951
16 -19800.9634 -27184.5711
17 -26826.5523 -19800.9634
18 -33396.5929 -26826.5523
19 -36576.7320 -33396.5929
20 -40780.3338 -36576.7320
21 -27939.4803 -40780.3338
22 -27030.5232 -27939.4803
23 -36841.9030 -27030.5232
24 -16651.7724 -36841.9030
25 -11675.7754 -16651.7724
26 -28175.6067 -11675.7754
27 -13123.9768 -28175.6067
28 -20057.7460 -13123.9768
29 -19650.8082 -20057.7460
30 -6754.6692 -19650.8082
31 -20104.6180 -6754.6692
32 -10366.7517 -20104.6180
33 -11548.1364 -10366.7517
34 -9395.6377 -11548.1364
35 -11002.6229 -9395.6377
36 13029.4034 -11002.6229
37 -4100.2421 13029.4034
38 -13116.4662 -4100.2421
39 2537.9464 -13116.4662
40 -6721.7186 2537.9464
41 -4809.8914 -6721.7186
42 11622.3662 -4809.8914
43 -2437.4141 11622.3662
44 1545.1185 -2437.4141
45 -6579.5156 1545.1185
46 -7666.2635 -6579.5156
47 110.3985 -7666.2635
48 6353.4213 110.3985
49 15633.2097 6353.4213
50 31771.2973 15633.2097
51 32946.4037 31771.2973
52 2415.5681 32946.4037
53 21220.1801 2415.5681
54 21938.5109 21220.1801
55 29548.2801 21938.5109
56 31554.8206 29548.2801
57 33916.5560 31554.8206
58 29186.8360 33916.5560
59 30197.9131 29186.8360
60 33650.5445 30197.9131
61 30980.6574 33650.5445
62 22255.2044 30980.6574
63 21538.0456 22255.2044
64 24852.2554 21538.0456
65 34155.6678 24852.2554
66 24518.6604 34155.6678
67 28905.4662 24518.6604
68 29074.7094 28905.4662
69 18397.0810 29074.7094
70 21413.6429 18397.0810
71 35463.5748 21413.6429
72 16047.6167 35463.5748
73 19121.1935 16047.6167
74 30259.5197 19121.1935
75 13793.7664 30259.5197
76 28927.8286 13793.7664
77 30158.9112 28927.8286
78 22269.7380 30158.9112
79 16528.4755 22269.7380
80 22536.2427 16528.4755
81 21168.5461 22536.2427
82 10119.4902 21168.5461
83 11905.0480 10119.4902
84 NA 11905.0480
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -18576.3580 -39106.5796
[2,] -20820.6534 -18576.3580
[3,] -30507.6143 -20820.6534
[4,] -9615.2241 -30507.6143
[5,] -34247.5071 -9615.2241
[6,] -40198.0135 -34247.5071
[7,] -15863.4576 -40198.0135
[8,] -33563.8057 -15863.4576
[9,] -27415.0509 -33563.8057
[10,] -16627.5448 -27415.0509
[11,] -29832.4085 -16627.5448
[12,] -13322.6338 -29832.4085
[13,] -31382.6851 -13322.6338
[14,] -22173.2951 -31382.6851
[15,] -27184.5711 -22173.2951
[16,] -19800.9634 -27184.5711
[17,] -26826.5523 -19800.9634
[18,] -33396.5929 -26826.5523
[19,] -36576.7320 -33396.5929
[20,] -40780.3338 -36576.7320
[21,] -27939.4803 -40780.3338
[22,] -27030.5232 -27939.4803
[23,] -36841.9030 -27030.5232
[24,] -16651.7724 -36841.9030
[25,] -11675.7754 -16651.7724
[26,] -28175.6067 -11675.7754
[27,] -13123.9768 -28175.6067
[28,] -20057.7460 -13123.9768
[29,] -19650.8082 -20057.7460
[30,] -6754.6692 -19650.8082
[31,] -20104.6180 -6754.6692
[32,] -10366.7517 -20104.6180
[33,] -11548.1364 -10366.7517
[34,] -9395.6377 -11548.1364
[35,] -11002.6229 -9395.6377
[36,] 13029.4034 -11002.6229
[37,] -4100.2421 13029.4034
[38,] -13116.4662 -4100.2421
[39,] 2537.9464 -13116.4662
[40,] -6721.7186 2537.9464
[41,] -4809.8914 -6721.7186
[42,] 11622.3662 -4809.8914
[43,] -2437.4141 11622.3662
[44,] 1545.1185 -2437.4141
[45,] -6579.5156 1545.1185
[46,] -7666.2635 -6579.5156
[47,] 110.3985 -7666.2635
[48,] 6353.4213 110.3985
[49,] 15633.2097 6353.4213
[50,] 31771.2973 15633.2097
[51,] 32946.4037 31771.2973
[52,] 2415.5681 32946.4037
[53,] 21220.1801 2415.5681
[54,] 21938.5109 21220.1801
[55,] 29548.2801 21938.5109
[56,] 31554.8206 29548.2801
[57,] 33916.5560 31554.8206
[58,] 29186.8360 33916.5560
[59,] 30197.9131 29186.8360
[60,] 33650.5445 30197.9131
[61,] 30980.6574 33650.5445
[62,] 22255.2044 30980.6574
[63,] 21538.0456 22255.2044
[64,] 24852.2554 21538.0456
[65,] 34155.6678 24852.2554
[66,] 24518.6604 34155.6678
[67,] 28905.4662 24518.6604
[68,] 29074.7094 28905.4662
[69,] 18397.0810 29074.7094
[70,] 21413.6429 18397.0810
[71,] 35463.5748 21413.6429
[72,] 16047.6167 35463.5748
[73,] 19121.1935 16047.6167
[74,] 30259.5197 19121.1935
[75,] 13793.7664 30259.5197
[76,] 28927.8286 13793.7664
[77,] 30158.9112 28927.8286
[78,] 22269.7380 30158.9112
[79,] 16528.4755 22269.7380
[80,] 22536.2427 16528.4755
[81,] 21168.5461 22536.2427
[82,] 10119.4902 21168.5461
[83,] 11905.0480 10119.4902
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -18576.3580 -39106.5796
2 -20820.6534 -18576.3580
3 -30507.6143 -20820.6534
4 -9615.2241 -30507.6143
5 -34247.5071 -9615.2241
6 -40198.0135 -34247.5071
7 -15863.4576 -40198.0135
8 -33563.8057 -15863.4576
9 -27415.0509 -33563.8057
10 -16627.5448 -27415.0509
11 -29832.4085 -16627.5448
12 -13322.6338 -29832.4085
13 -31382.6851 -13322.6338
14 -22173.2951 -31382.6851
15 -27184.5711 -22173.2951
16 -19800.9634 -27184.5711
17 -26826.5523 -19800.9634
18 -33396.5929 -26826.5523
19 -36576.7320 -33396.5929
20 -40780.3338 -36576.7320
21 -27939.4803 -40780.3338
22 -27030.5232 -27939.4803
23 -36841.9030 -27030.5232
24 -16651.7724 -36841.9030
25 -11675.7754 -16651.7724
26 -28175.6067 -11675.7754
27 -13123.9768 -28175.6067
28 -20057.7460 -13123.9768
29 -19650.8082 -20057.7460
30 -6754.6692 -19650.8082
31 -20104.6180 -6754.6692
32 -10366.7517 -20104.6180
33 -11548.1364 -10366.7517
34 -9395.6377 -11548.1364
35 -11002.6229 -9395.6377
36 13029.4034 -11002.6229
37 -4100.2421 13029.4034
38 -13116.4662 -4100.2421
39 2537.9464 -13116.4662
40 -6721.7186 2537.9464
41 -4809.8914 -6721.7186
42 11622.3662 -4809.8914
43 -2437.4141 11622.3662
44 1545.1185 -2437.4141
45 -6579.5156 1545.1185
46 -7666.2635 -6579.5156
47 110.3985 -7666.2635
48 6353.4213 110.3985
49 15633.2097 6353.4213
50 31771.2973 15633.2097
51 32946.4037 31771.2973
52 2415.5681 32946.4037
53 21220.1801 2415.5681
54 21938.5109 21220.1801
55 29548.2801 21938.5109
56 31554.8206 29548.2801
57 33916.5560 31554.8206
58 29186.8360 33916.5560
59 30197.9131 29186.8360
60 33650.5445 30197.9131
61 30980.6574 33650.5445
62 22255.2044 30980.6574
63 21538.0456 22255.2044
64 24852.2554 21538.0456
65 34155.6678 24852.2554
66 24518.6604 34155.6678
67 28905.4662 24518.6604
68 29074.7094 28905.4662
69 18397.0810 29074.7094
70 21413.6429 18397.0810
71 35463.5748 21413.6429
72 16047.6167 35463.5748
73 19121.1935 16047.6167
74 30259.5197 19121.1935
75 13793.7664 30259.5197
76 28927.8286 13793.7664
77 30158.9112 28927.8286
78 22269.7380 30158.9112
79 16528.4755 22269.7380
80 22536.2427 16528.4755
81 21168.5461 22536.2427
82 10119.4902 21168.5461
83 11905.0480 10119.4902
> 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/76rh01229727827.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/89frj1229727827.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/9s56d1229727828.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/10rb8t1229727828.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/11m6lz1229727828.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/12onw01229727828.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/13ewjg1229727828.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/14wlnb1229727828.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/155aen1229727828.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/16ds061229727828.tab")
+ }
>
> system("convert tmp/1rcg51229727827.ps tmp/1rcg51229727827.png")
> system("convert tmp/2619j1229727827.ps tmp/2619j1229727827.png")
> system("convert tmp/39fg81229727827.ps tmp/39fg81229727827.png")
> system("convert tmp/43quo1229727827.ps tmp/43quo1229727827.png")
> system("convert tmp/57h281229727827.ps tmp/57h281229727827.png")
> system("convert tmp/6xhev1229727827.ps tmp/6xhev1229727827.png")
> system("convert tmp/76rh01229727827.ps tmp/76rh01229727827.png")
> system("convert tmp/89frj1229727827.ps tmp/89frj1229727827.png")
> system("convert tmp/9s56d1229727828.ps tmp/9s56d1229727828.png")
> system("convert tmp/10rb8t1229727828.ps tmp/10rb8t1229727828.png")
>
>
> proc.time()
user system elapsed
2.986 1.705 3.842