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
+ ,11554.5
+ ,173666
+ ,13182.1
+ ,165688
+ ,14800.1
+ ,161570
+ ,12150.7
+ ,156145
+ ,14478.2
+ ,153730
+ ,13253.9
+ ,182698
+ ,12036.8
+ ,200765
+ ,12653.2
+ ,176512
+ ,14035.4
+ ,166618
+ ,14571.4
+ ,158644
+ ,15400.9
+ ,159585
+ ,14283.2
+ ,163095
+ ,14485.3
+ ,159044
+ ,14196.3
+ ,155511
+ ,15559.1
+ ,153745
+ ,13767.4
+ ,150569
+ ,14634
+ ,150605
+ ,14381.1
+ ,179612
+ ,12509.9
+ ,194690
+ ,12122.3
+ ,189917
+ ,13122.3
+ ,184128
+ ,13908.7
+ ,175335
+ ,13456.5
+ ,179566
+ ,12441.6
+ ,181140
+ ,12953
+ ,177876
+ ,13057.2
+ ,175041
+ ,14350.1
+ ,169292
+ ,13830.2
+ ,166070
+ ,13755.5
+ ,166972
+ ,13574.4
+ ,206348
+ ,12802.6
+ ,215706
+ ,11737.3
+ ,202108
+ ,13850.2
+ ,195411
+ ,15081.8
+ ,193111
+ ,13653.3
+ ,195198
+ ,14019.1
+ ,198770
+ ,13962
+ ,194163
+ ,13768.7
+ ,190420
+ ,14747.1
+ ,189733
+ ,13858.1
+ ,186029
+ ,13188
+ ,191531
+ ,13693.1
+ ,232571
+ ,12970
+ ,243477
+ ,11392.8
+ ,227247
+ ,13985.2
+ ,217859
+ ,14994.7
+ ,208679
+ ,13584.7
+ ,213188
+ ,14257.8
+ ,216234
+ ,13553.4
+ ,213586
+ ,14007.3
+ ,209465
+ ,16535.8
+ ,204045
+ ,14721.4
+ ,200237
+ ,13664.6
+ ,203666
+ ,16405.9
+ ,241476
+ ,13829.4
+ ,260307
+ ,13735.6
+ ,243324
+ ,15870.5
+ ,244460
+ ,15962.4
+ ,233575
+ ,15744.1
+ ,237217
+ ,16083.7
+ ,235243
+ ,14863.9
+ ,230354
+ ,15533.1
+ ,227184
+ ,17473.1
+ ,221678
+ ,15925.5
+ ,217142
+ ,15573.7
+ ,219452
+ ,17495
+ ,256446
+ ,14155.8
+ ,265845
+ ,14913.9
+ ,248624
+ ,17250.4
+ ,241114
+ ,15879.8
+ ,229245
+ ,17647.8
+ ,231805
+ ,17749.9
+ ,219277
+ ,17111.8
+ ,219313
+ ,16934.8
+ ,212610
+ ,20280
+ ,214771
+ ,16238.2
+ ,211142
+ ,17896.1
+ ,211457
+ ,18089.3
+ ,240048
+ ,15660
+ ,240636
+ ,16162.4
+ ,230580
+ ,17850.1
+ ,208795
+ ,18520.4
+ ,197922
+ ,18524.7
+ ,194596
+ ,16843.7)
+ ,dim=c(2
+ ,84)
+ ,dimnames=list(c('werkloosheid'
+ ,'invoer')
+ ,1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('werkloosheid','invoer'),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 invoer M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 180144 11554.5 1 0 0 0 0 0 0 0 0 0 0
2 173666 13182.1 0 1 0 0 0 0 0 0 0 0 0
3 165688 14800.1 0 0 1 0 0 0 0 0 0 0 0
4 161570 12150.7 0 0 0 1 0 0 0 0 0 0 0
5 156145 14478.2 0 0 0 0 1 0 0 0 0 0 0
6 153730 13253.9 0 0 0 0 0 1 0 0 0 0 0
7 182698 12036.8 0 0 0 0 0 0 1 0 0 0 0
8 200765 12653.2 0 0 0 0 0 0 0 1 0 0 0
9 176512 14035.4 0 0 0 0 0 0 0 0 1 0 0
10 166618 14571.4 0 0 0 0 0 0 0 0 0 1 0
11 158644 15400.9 0 0 0 0 0 0 0 0 0 0 1
12 159585 14283.2 0 0 0 0 0 0 0 0 0 0 0
13 163095 14485.3 1 0 0 0 0 0 0 0 0 0 0
14 159044 14196.3 0 1 0 0 0 0 0 0 0 0 0
15 155511 15559.1 0 0 1 0 0 0 0 0 0 0 0
16 153745 13767.4 0 0 0 1 0 0 0 0 0 0 0
17 150569 14634.0 0 0 0 0 1 0 0 0 0 0 0
18 150605 14381.1 0 0 0 0 0 1 0 0 0 0 0
19 179612 12509.9 0 0 0 0 0 0 1 0 0 0 0
20 194690 12122.3 0 0 0 0 0 0 0 1 0 0 0
21 189917 13122.3 0 0 0 0 0 0 0 0 1 0 0
22 184128 13908.7 0 0 0 0 0 0 0 0 0 1 0
23 175335 13456.5 0 0 0 0 0 0 0 0 0 0 1
24 179566 12441.6 0 0 0 0 0 0 0 0 0 0 0
25 181140 12953.0 1 0 0 0 0 0 0 0 0 0 0
26 177876 13057.2 0 1 0 0 0 0 0 0 0 0 0
27 175041 14350.1 0 0 1 0 0 0 0 0 0 0 0
28 169292 13830.2 0 0 0 1 0 0 0 0 0 0 0
29 166070 13755.5 0 0 0 0 1 0 0 0 0 0 0
30 166972 13574.4 0 0 0 0 0 1 0 0 0 0 0
31 206348 12802.6 0 0 0 0 0 0 1 0 0 0 0
32 215706 11737.3 0 0 0 0 0 0 0 1 0 0 0
33 202108 13850.2 0 0 0 0 0 0 0 0 1 0 0
34 195411 15081.8 0 0 0 0 0 0 0 0 0 1 0
35 193111 13653.3 0 0 0 0 0 0 0 0 0 0 1
36 195198 14019.1 0 0 0 0 0 0 0 0 0 0 0
37 198770 13962.0 1 0 0 0 0 0 0 0 0 0 0
38 194163 13768.7 0 1 0 0 0 0 0 0 0 0 0
39 190420 14747.1 0 0 1 0 0 0 0 0 0 0 0
40 189733 13858.1 0 0 0 1 0 0 0 0 0 0 0
41 186029 13188.0 0 0 0 0 1 0 0 0 0 0 0
42 191531 13693.1 0 0 0 0 0 1 0 0 0 0 0
43 232571 12970.0 0 0 0 0 0 0 1 0 0 0 0
44 243477 11392.8 0 0 0 0 0 0 0 1 0 0 0
45 227247 13985.2 0 0 0 0 0 0 0 0 1 0 0
46 217859 14994.7 0 0 0 0 0 0 0 0 0 1 0
47 208679 13584.7 0 0 0 0 0 0 0 0 0 0 1
48 213188 14257.8 0 0 0 0 0 0 0 0 0 0 0
49 216234 13553.4 1 0 0 0 0 0 0 0 0 0 0
50 213586 14007.3 0 1 0 0 0 0 0 0 0 0 0
51 209465 16535.8 0 0 1 0 0 0 0 0 0 0 0
52 204045 14721.4 0 0 0 1 0 0 0 0 0 0 0
53 200237 13664.6 0 0 0 0 1 0 0 0 0 0 0
54 203666 16405.9 0 0 0 0 0 1 0 0 0 0 0
55 241476 13829.4 0 0 0 0 0 0 1 0 0 0 0
56 260307 13735.6 0 0 0 0 0 0 0 1 0 0 0
57 243324 15870.5 0 0 0 0 0 0 0 0 1 0 0
58 244460 15962.4 0 0 0 0 0 0 0 0 0 1 0
59 233575 15744.1 0 0 0 0 0 0 0 0 0 0 1
60 237217 16083.7 0 0 0 0 0 0 0 0 0 0 0
61 235243 14863.9 1 0 0 0 0 0 0 0 0 0 0
62 230354 15533.1 0 1 0 0 0 0 0 0 0 0 0
63 227184 17473.1 0 0 1 0 0 0 0 0 0 0 0
64 221678 15925.5 0 0 0 1 0 0 0 0 0 0 0
65 217142 15573.7 0 0 0 0 1 0 0 0 0 0 0
66 219452 17495.0 0 0 0 0 0 1 0 0 0 0 0
67 256446 14155.8 0 0 0 0 0 0 1 0 0 0 0
68 265845 14913.9 0 0 0 0 0 0 0 1 0 0 0
69 248624 17250.4 0 0 0 0 0 0 0 0 1 0 0
70 241114 15879.8 0 0 0 0 0 0 0 0 0 1 0
71 229245 17647.8 0 0 0 0 0 0 0 0 0 0 1
72 231805 17749.9 0 0 0 0 0 0 0 0 0 0 0
73 219277 17111.8 1 0 0 0 0 0 0 0 0 0 0
74 219313 16934.8 0 1 0 0 0 0 0 0 0 0 0
75 212610 20280.0 0 0 1 0 0 0 0 0 0 0 0
76 214771 16238.2 0 0 0 1 0 0 0 0 0 0 0
77 211142 17896.1 0 0 0 0 1 0 0 0 0 0 0
78 211457 18089.3 0 0 0 0 0 1 0 0 0 0 0
79 240048 15660.0 0 0 0 0 0 0 1 0 0 0 0
80 240636 16162.4 0 0 0 0 0 0 0 1 0 0 0
81 230580 17850.1 0 0 0 0 0 0 0 0 1 0 0
82 208795 18520.4 0 0 0 0 0 0 0 0 0 1 0
83 197922 18524.7 0 0 0 0 0 0 0 0 0 0 1
84 194596 16843.7 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) invoer M1 M2 M3 M4
54343.510 9.754 7560.873 801.444 -21987.349 -6532.030
M5 M6 M7 M8 M9 M10
-14220.753 -17939.992 34614.645 48098.893 14910.893 2232.342
M11
-5342.735
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-40571.0 -12936.2 -160.6 17708.2 32193.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 54343.510 24161.850 2.249 0.027606 *
invoer 9.754 1.502 6.492 9.9e-09 ***
M1 7560.873 11871.646 0.637 0.526248
M2 801.444 11819.579 0.068 0.946130
M3 -21987.349 11897.416 -1.848 0.068756 .
M4 -6532.030 11823.311 -0.552 0.582361
M5 -14220.753 11782.882 -1.207 0.231477
M6 -17939.992 11773.648 -1.524 0.132015
M7 34614.645 12036.327 2.876 0.005316 **
M8 48098.893 12095.074 3.977 0.000166 ***
M9 14910.893 11770.924 1.267 0.209382
M10 2232.342 11791.294 0.189 0.850381
M11 -5342.735 11781.412 -0.453 0.651580
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 22020 on 71 degrees of freedom
Multiple R-squared: 0.5218, Adjusted R-squared: 0.4409
F-statistic: 6.455 on 12 and 71 DF, p-value: 1.311e-07
> 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,] 1.349044e-02 2.698088e-02 0.986509558
[2,] 3.901117e-03 7.802234e-03 0.996098883
[3,] 1.110162e-03 2.220323e-03 0.998889838
[4,] 2.661518e-04 5.323036e-04 0.999733848
[5,] 3.119415e-04 6.238831e-04 0.999688058
[6,] 1.770125e-04 3.540249e-04 0.999822988
[7,] 3.187236e-04 6.374472e-04 0.999681276
[8,] 1.068717e-04 2.137435e-04 0.999893128
[9,] 4.363026e-05 8.726053e-05 0.999956370
[10,] 4.081311e-05 8.162623e-05 0.999959187
[11,] 2.409826e-05 4.819653e-05 0.999975902
[12,] 1.570112e-05 3.140224e-05 0.999984299
[13,] 1.642253e-04 3.284507e-04 0.999835775
[14,] 1.301553e-04 2.603106e-04 0.999869845
[15,] 2.126771e-04 4.253542e-04 0.999787323
[16,] 5.914209e-03 1.182842e-02 0.994085791
[17,] 7.785142e-03 1.557028e-02 0.992214858
[18,] 1.898883e-02 3.797766e-02 0.981011171
[19,] 6.006613e-02 1.201323e-01 0.939933874
[20,] 8.308052e-02 1.661610e-01 0.916919483
[21,] 1.713784e-01 3.427567e-01 0.828621644
[22,] 3.027043e-01 6.054086e-01 0.697295697
[23,] 4.111918e-01 8.223837e-01 0.588808175
[24,] 4.910364e-01 9.820727e-01 0.508963636
[25,] 6.274388e-01 7.451225e-01 0.372561234
[26,] 6.894390e-01 6.211220e-01 0.310560979
[27,] 7.881776e-01 4.236448e-01 0.211822419
[28,] 8.905404e-01 2.189193e-01 0.109459629
[29,] 9.176982e-01 1.646036e-01 0.082301806
[30,] 9.441881e-01 1.116239e-01 0.055811938
[31,] 9.591037e-01 8.179263e-02 0.040896314
[32,] 9.644197e-01 7.116064e-02 0.035580322
[33,] 9.704955e-01 5.900904e-02 0.029504520
[34,] 9.765927e-01 4.681464e-02 0.023407321
[35,] 9.823714e-01 3.525725e-02 0.017628626
[36,] 9.871441e-01 2.571190e-02 0.012855949
[37,] 9.897590e-01 2.048203e-02 0.010241015
[38,] 9.966378e-01 6.724381e-03 0.003362191
[39,] 9.976404e-01 4.719224e-03 0.002359612
[40,] 9.976229e-01 4.754223e-03 0.002377111
[41,] 9.965937e-01 6.812512e-03 0.003406256
[42,] 9.950361e-01 9.927765e-03 0.004963883
[43,] 9.933624e-01 1.327520e-02 0.006637601
[44,] 9.884197e-01 2.316068e-02 0.011580338
[45,] 9.843481e-01 3.130376e-02 0.015651880
[46,] 9.721470e-01 5.570608e-02 0.027853041
[47,] 9.499725e-01 1.000549e-01 0.050027455
[48,] 9.174689e-01 1.650622e-01 0.082531120
[49,] 8.609933e-01 2.780135e-01 0.139006730
[50,] 8.160474e-01 3.679051e-01 0.183952561
[51,] 7.048555e-01 5.902890e-01 0.295144515
[52,] 5.702878e-01 8.594244e-01 0.429712218
[53,] 4.201477e-01 8.402954e-01 0.579852308
> postscript(file="/var/www/html/rcomp/tmp/1g30c1229718036.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/20gwp1229718036.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/3dmrg1229718036.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/43zbp1229718036.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/5y3ut1229718036.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 7
5541.675 -10051.852 -11022.379 -4754.517 -25192.292 -11946.724 -23662.258
8 9 10 11 12 13 14
-25091.623 -29638.045 -32081.422 -40570.955 -34071.094 -40093.169 -34565.951
15 16 17 18 19 20 21
-28602.360 -28348.158 -32287.903 -26065.979 -31362.685 -25988.438 -7327.035
22 23 24 25 26 27 28
-8107.713 -4915.060 3872.131 -7102.731 -4623.628 2719.740 -13413.684
29 30 31 32 33 34 35
-8218.367 -1830.752 -7481.563 -1217.303 -2235.679 -8266.659 10941.432
36 37 38 39 40 41 42
4117.831 685.889 4723.688 14226.561 6755.191 17275.799 21570.496
43 44 45 46 47 48 49
17108.685 29913.811 21586.586 15030.880 27178.529 19779.648 22135.209
50 51 52 53 54 55 56
21819.479 15825.302 12646.910 26835.235 7245.936 17631.443 23893.083
57 58 59 60 61 62 63
19275.128 32193.323 31012.610 25999.554 28362.119 23705.440 24402.255
64 65 66 67 68 69 70
18535.603 25119.642 12409.293 29417.869 17938.419 11116.139 29652.971
71 72 73 74 75 76 77
8114.687 4336.110 -9528.992 -1007.177 -17549.118 8578.654 -3532.113
78 79 80 81 82 83 84
-1382.270 -1651.492 -19447.948 -12777.093 -28421.379 -31761.243 -24034.180
> postscript(file="/var/www/html/rcomp/tmp/6vikm1229718036.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 5541.675 NA
1 -10051.852 5541.675
2 -11022.379 -10051.852
3 -4754.517 -11022.379
4 -25192.292 -4754.517
5 -11946.724 -25192.292
6 -23662.258 -11946.724
7 -25091.623 -23662.258
8 -29638.045 -25091.623
9 -32081.422 -29638.045
10 -40570.955 -32081.422
11 -34071.094 -40570.955
12 -40093.169 -34071.094
13 -34565.951 -40093.169
14 -28602.360 -34565.951
15 -28348.158 -28602.360
16 -32287.903 -28348.158
17 -26065.979 -32287.903
18 -31362.685 -26065.979
19 -25988.438 -31362.685
20 -7327.035 -25988.438
21 -8107.713 -7327.035
22 -4915.060 -8107.713
23 3872.131 -4915.060
24 -7102.731 3872.131
25 -4623.628 -7102.731
26 2719.740 -4623.628
27 -13413.684 2719.740
28 -8218.367 -13413.684
29 -1830.752 -8218.367
30 -7481.563 -1830.752
31 -1217.303 -7481.563
32 -2235.679 -1217.303
33 -8266.659 -2235.679
34 10941.432 -8266.659
35 4117.831 10941.432
36 685.889 4117.831
37 4723.688 685.889
38 14226.561 4723.688
39 6755.191 14226.561
40 17275.799 6755.191
41 21570.496 17275.799
42 17108.685 21570.496
43 29913.811 17108.685
44 21586.586 29913.811
45 15030.880 21586.586
46 27178.529 15030.880
47 19779.648 27178.529
48 22135.209 19779.648
49 21819.479 22135.209
50 15825.302 21819.479
51 12646.910 15825.302
52 26835.235 12646.910
53 7245.936 26835.235
54 17631.443 7245.936
55 23893.083 17631.443
56 19275.128 23893.083
57 32193.323 19275.128
58 31012.610 32193.323
59 25999.554 31012.610
60 28362.119 25999.554
61 23705.440 28362.119
62 24402.255 23705.440
63 18535.603 24402.255
64 25119.642 18535.603
65 12409.293 25119.642
66 29417.869 12409.293
67 17938.419 29417.869
68 11116.139 17938.419
69 29652.971 11116.139
70 8114.687 29652.971
71 4336.110 8114.687
72 -9528.992 4336.110
73 -1007.177 -9528.992
74 -17549.118 -1007.177
75 8578.654 -17549.118
76 -3532.113 8578.654
77 -1382.270 -3532.113
78 -1651.492 -1382.270
79 -19447.948 -1651.492
80 -12777.093 -19447.948
81 -28421.379 -12777.093
82 -31761.243 -28421.379
83 -24034.180 -31761.243
84 NA -24034.180
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -10051.852 5541.675
[2,] -11022.379 -10051.852
[3,] -4754.517 -11022.379
[4,] -25192.292 -4754.517
[5,] -11946.724 -25192.292
[6,] -23662.258 -11946.724
[7,] -25091.623 -23662.258
[8,] -29638.045 -25091.623
[9,] -32081.422 -29638.045
[10,] -40570.955 -32081.422
[11,] -34071.094 -40570.955
[12,] -40093.169 -34071.094
[13,] -34565.951 -40093.169
[14,] -28602.360 -34565.951
[15,] -28348.158 -28602.360
[16,] -32287.903 -28348.158
[17,] -26065.979 -32287.903
[18,] -31362.685 -26065.979
[19,] -25988.438 -31362.685
[20,] -7327.035 -25988.438
[21,] -8107.713 -7327.035
[22,] -4915.060 -8107.713
[23,] 3872.131 -4915.060
[24,] -7102.731 3872.131
[25,] -4623.628 -7102.731
[26,] 2719.740 -4623.628
[27,] -13413.684 2719.740
[28,] -8218.367 -13413.684
[29,] -1830.752 -8218.367
[30,] -7481.563 -1830.752
[31,] -1217.303 -7481.563
[32,] -2235.679 -1217.303
[33,] -8266.659 -2235.679
[34,] 10941.432 -8266.659
[35,] 4117.831 10941.432
[36,] 685.889 4117.831
[37,] 4723.688 685.889
[38,] 14226.561 4723.688
[39,] 6755.191 14226.561
[40,] 17275.799 6755.191
[41,] 21570.496 17275.799
[42,] 17108.685 21570.496
[43,] 29913.811 17108.685
[44,] 21586.586 29913.811
[45,] 15030.880 21586.586
[46,] 27178.529 15030.880
[47,] 19779.648 27178.529
[48,] 22135.209 19779.648
[49,] 21819.479 22135.209
[50,] 15825.302 21819.479
[51,] 12646.910 15825.302
[52,] 26835.235 12646.910
[53,] 7245.936 26835.235
[54,] 17631.443 7245.936
[55,] 23893.083 17631.443
[56,] 19275.128 23893.083
[57,] 32193.323 19275.128
[58,] 31012.610 32193.323
[59,] 25999.554 31012.610
[60,] 28362.119 25999.554
[61,] 23705.440 28362.119
[62,] 24402.255 23705.440
[63,] 18535.603 24402.255
[64,] 25119.642 18535.603
[65,] 12409.293 25119.642
[66,] 29417.869 12409.293
[67,] 17938.419 29417.869
[68,] 11116.139 17938.419
[69,] 29652.971 11116.139
[70,] 8114.687 29652.971
[71,] 4336.110 8114.687
[72,] -9528.992 4336.110
[73,] -1007.177 -9528.992
[74,] -17549.118 -1007.177
[75,] 8578.654 -17549.118
[76,] -3532.113 8578.654
[77,] -1382.270 -3532.113
[78,] -1651.492 -1382.270
[79,] -19447.948 -1651.492
[80,] -12777.093 -19447.948
[81,] -28421.379 -12777.093
[82,] -31761.243 -28421.379
[83,] -24034.180 -31761.243
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -10051.852 5541.675
2 -11022.379 -10051.852
3 -4754.517 -11022.379
4 -25192.292 -4754.517
5 -11946.724 -25192.292
6 -23662.258 -11946.724
7 -25091.623 -23662.258
8 -29638.045 -25091.623
9 -32081.422 -29638.045
10 -40570.955 -32081.422
11 -34071.094 -40570.955
12 -40093.169 -34071.094
13 -34565.951 -40093.169
14 -28602.360 -34565.951
15 -28348.158 -28602.360
16 -32287.903 -28348.158
17 -26065.979 -32287.903
18 -31362.685 -26065.979
19 -25988.438 -31362.685
20 -7327.035 -25988.438
21 -8107.713 -7327.035
22 -4915.060 -8107.713
23 3872.131 -4915.060
24 -7102.731 3872.131
25 -4623.628 -7102.731
26 2719.740 -4623.628
27 -13413.684 2719.740
28 -8218.367 -13413.684
29 -1830.752 -8218.367
30 -7481.563 -1830.752
31 -1217.303 -7481.563
32 -2235.679 -1217.303
33 -8266.659 -2235.679
34 10941.432 -8266.659
35 4117.831 10941.432
36 685.889 4117.831
37 4723.688 685.889
38 14226.561 4723.688
39 6755.191 14226.561
40 17275.799 6755.191
41 21570.496 17275.799
42 17108.685 21570.496
43 29913.811 17108.685
44 21586.586 29913.811
45 15030.880 21586.586
46 27178.529 15030.880
47 19779.648 27178.529
48 22135.209 19779.648
49 21819.479 22135.209
50 15825.302 21819.479
51 12646.910 15825.302
52 26835.235 12646.910
53 7245.936 26835.235
54 17631.443 7245.936
55 23893.083 17631.443
56 19275.128 23893.083
57 32193.323 19275.128
58 31012.610 32193.323
59 25999.554 31012.610
60 28362.119 25999.554
61 23705.440 28362.119
62 24402.255 23705.440
63 18535.603 24402.255
64 25119.642 18535.603
65 12409.293 25119.642
66 29417.869 12409.293
67 17938.419 29417.869
68 11116.139 17938.419
69 29652.971 11116.139
70 8114.687 29652.971
71 4336.110 8114.687
72 -9528.992 4336.110
73 -1007.177 -9528.992
74 -17549.118 -1007.177
75 8578.654 -17549.118
76 -3532.113 8578.654
77 -1382.270 -3532.113
78 -1651.492 -1382.270
79 -19447.948 -1651.492
80 -12777.093 -19447.948
81 -28421.379 -12777.093
82 -31761.243 -28421.379
83 -24034.180 -31761.243
> 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/72ihi1229718036.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/8jud81229718036.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/9b8dr1229718036.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/10xxf71229718036.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/11pykd1229718036.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/12y0g61229718036.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/13p72t1229718036.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/14pds81229718036.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/1595yc1229718036.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/16hhgx1229718036.tab")
+ }
>
> system("convert tmp/1g30c1229718036.ps tmp/1g30c1229718036.png")
> system("convert tmp/20gwp1229718036.ps tmp/20gwp1229718036.png")
> system("convert tmp/3dmrg1229718036.ps tmp/3dmrg1229718036.png")
> system("convert tmp/43zbp1229718036.ps tmp/43zbp1229718036.png")
> system("convert tmp/5y3ut1229718036.ps tmp/5y3ut1229718036.png")
> system("convert tmp/6vikm1229718036.ps tmp/6vikm1229718036.png")
> system("convert tmp/72ihi1229718036.ps tmp/72ihi1229718036.png")
> system("convert tmp/8jud81229718036.ps tmp/8jud81229718036.png")
> system("convert tmp/9b8dr1229718036.ps tmp/9b8dr1229718036.png")
> system("convert tmp/10xxf71229718036.ps tmp/10xxf71229718036.png")
>
>
> proc.time()
user system elapsed
2.809 1.658 3.609