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(11554.5
+ ,180144
+ ,13182.1
+ ,173666
+ ,14800.1
+ ,165688
+ ,12150.7
+ ,161570
+ ,14478.2
+ ,156145
+ ,13253.9
+ ,153730
+ ,12036.8
+ ,182698
+ ,12653.2
+ ,200765
+ ,14035.4
+ ,176512
+ ,14571.4
+ ,166618
+ ,15400.9
+ ,158644
+ ,14283.2
+ ,159585
+ ,14485.3
+ ,163095
+ ,14196.3
+ ,159044
+ ,15559.1
+ ,155511
+ ,13767.4
+ ,153745
+ ,14634
+ ,150569
+ ,14381.1
+ ,150605
+ ,12509.9
+ ,179612
+ ,12122.3
+ ,194690
+ ,13122.3
+ ,189917
+ ,13908.7
+ ,184128
+ ,13456.5
+ ,175335
+ ,12441.6
+ ,179566
+ ,12953
+ ,181140
+ ,13057.2
+ ,177876
+ ,14350.1
+ ,175041
+ ,13830.2
+ ,169292
+ ,13755.5
+ ,166070
+ ,13574.4
+ ,166972
+ ,12802.6
+ ,206348
+ ,11737.3
+ ,215706
+ ,13850.2
+ ,202108
+ ,15081.8
+ ,195411
+ ,13653.3
+ ,193111
+ ,14019.1
+ ,195198
+ ,13962
+ ,198770
+ ,13768.7
+ ,194163
+ ,14747.1
+ ,190420
+ ,13858.1
+ ,189733
+ ,13188
+ ,186029
+ ,13693.1
+ ,191531
+ ,12970
+ ,232571
+ ,11392.8
+ ,243477
+ ,13985.2
+ ,227247
+ ,14994.7
+ ,217859
+ ,13584.7
+ ,208679
+ ,14257.8
+ ,213188
+ ,13553.4
+ ,216234
+ ,14007.3
+ ,213586
+ ,16535.8
+ ,209465
+ ,14721.4
+ ,204045
+ ,13664.6
+ ,200237
+ ,16805.9
+ ,203666
+ ,13829.4
+ ,241476
+ ,13735.6
+ ,260307
+ ,15870.5
+ ,243324
+ ,15962.4
+ ,244460
+ ,15744.1
+ ,233575
+ ,16083.7
+ ,237217
+ ,14863.9
+ ,235243
+ ,15533.1
+ ,230354
+ ,17473.1
+ ,227184
+ ,15925.5
+ ,221678
+ ,15573.7
+ ,217142
+ ,17495
+ ,219452
+ ,14155.8
+ ,256446
+ ,14913.9
+ ,265845
+ ,17250.4
+ ,248624
+ ,15879.8
+ ,241114
+ ,17647.8
+ ,229245
+ ,17749.9
+ ,231805
+ ,17111.8
+ ,219277
+ ,16934.8
+ ,219313
+ ,20280
+ ,212610
+ ,16238.2
+ ,214771
+ ,17896.1
+ ,211142
+ ,18089.3
+ ,211457
+ ,15660
+ ,240048
+ ,16162.4
+ ,240636
+ ,17850.1
+ ,230580
+ ,18520.4
+ ,208795
+ ,18524.7
+ ,197922
+ ,16843.7
+ ,194596)
+ ,dim=c(2
+ ,84)
+ ,dimnames=list(c('invoer'
+ ,'werkloosheid')
+ ,1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('invoer','werkloosheid'),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
invoer werkloosheid M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 11554.5 180144 1 0 0 0 0 0 0 0 0 0 0
2 13182.1 173666 0 1 0 0 0 0 0 0 0 0 0
3 14800.1 165688 0 0 1 0 0 0 0 0 0 0 0
4 12150.7 161570 0 0 0 1 0 0 0 0 0 0 0
5 14478.2 156145 0 0 0 0 1 0 0 0 0 0 0
6 13253.9 153730 0 0 0 0 0 1 0 0 0 0 0
7 12036.8 182698 0 0 0 0 0 0 1 0 0 0 0
8 12653.2 200765 0 0 0 0 0 0 0 1 0 0 0
9 14035.4 176512 0 0 0 0 0 0 0 0 1 0 0
10 14571.4 166618 0 0 0 0 0 0 0 0 0 1 0
11 15400.9 158644 0 0 0 0 0 0 0 0 0 0 1
12 14283.2 159585 0 0 0 0 0 0 0 0 0 0 0
13 14485.3 163095 1 0 0 0 0 0 0 0 0 0 0
14 14196.3 159044 0 1 0 0 0 0 0 0 0 0 0
15 15559.1 155511 0 0 1 0 0 0 0 0 0 0 0
16 13767.4 153745 0 0 0 1 0 0 0 0 0 0 0
17 14634.0 150569 0 0 0 0 1 0 0 0 0 0 0
18 14381.1 150605 0 0 0 0 0 1 0 0 0 0 0
19 12509.9 179612 0 0 0 0 0 0 1 0 0 0 0
20 12122.3 194690 0 0 0 0 0 0 0 1 0 0 0
21 13122.3 189917 0 0 0 0 0 0 0 0 1 0 0
22 13908.7 184128 0 0 0 0 0 0 0 0 0 1 0
23 13456.5 175335 0 0 0 0 0 0 0 0 0 0 1
24 12441.6 179566 0 0 0 0 0 0 0 0 0 0 0
25 12953.0 181140 1 0 0 0 0 0 0 0 0 0 0
26 13057.2 177876 0 1 0 0 0 0 0 0 0 0 0
27 14350.1 175041 0 0 1 0 0 0 0 0 0 0 0
28 13830.2 169292 0 0 0 1 0 0 0 0 0 0 0
29 13755.5 166070 0 0 0 0 1 0 0 0 0 0 0
30 13574.4 166972 0 0 0 0 0 1 0 0 0 0 0
31 12802.6 206348 0 0 0 0 0 0 1 0 0 0 0
32 11737.3 215706 0 0 0 0 0 0 0 1 0 0 0
33 13850.2 202108 0 0 0 0 0 0 0 0 1 0 0
34 15081.8 195411 0 0 0 0 0 0 0 0 0 1 0
35 13653.3 193111 0 0 0 0 0 0 0 0 0 0 1
36 14019.1 195198 0 0 0 0 0 0 0 0 0 0 0
37 13962.0 198770 1 0 0 0 0 0 0 0 0 0 0
38 13768.7 194163 0 1 0 0 0 0 0 0 0 0 0
39 14747.1 190420 0 0 1 0 0 0 0 0 0 0 0
40 13858.1 189733 0 0 0 1 0 0 0 0 0 0 0
41 13188.0 186029 0 0 0 0 1 0 0 0 0 0 0
42 13693.1 191531 0 0 0 0 0 1 0 0 0 0 0
43 12970.0 232571 0 0 0 0 0 0 1 0 0 0 0
44 11392.8 243477 0 0 0 0 0 0 0 1 0 0 0
45 13985.2 227247 0 0 0 0 0 0 0 0 1 0 0
46 14994.7 217859 0 0 0 0 0 0 0 0 0 1 0
47 13584.7 208679 0 0 0 0 0 0 0 0 0 0 1
48 14257.8 213188 0 0 0 0 0 0 0 0 0 0 0
49 13553.4 216234 1 0 0 0 0 0 0 0 0 0 0
50 14007.3 213586 0 1 0 0 0 0 0 0 0 0 0
51 16535.8 209465 0 0 1 0 0 0 0 0 0 0 0
52 14721.4 204045 0 0 0 1 0 0 0 0 0 0 0
53 13664.6 200237 0 0 0 0 1 0 0 0 0 0 0
54 16805.9 203666 0 0 0 0 0 1 0 0 0 0 0
55 13829.4 241476 0 0 0 0 0 0 1 0 0 0 0
56 13735.6 260307 0 0 0 0 0 0 0 1 0 0 0
57 15870.5 243324 0 0 0 0 0 0 0 0 1 0 0
58 15962.4 244460 0 0 0 0 0 0 0 0 0 1 0
59 15744.1 233575 0 0 0 0 0 0 0 0 0 0 1
60 16083.7 237217 0 0 0 0 0 0 0 0 0 0 0
61 14863.9 235243 1 0 0 0 0 0 0 0 0 0 0
62 15533.1 230354 0 1 0 0 0 0 0 0 0 0 0
63 17473.1 227184 0 0 1 0 0 0 0 0 0 0 0
64 15925.5 221678 0 0 0 1 0 0 0 0 0 0 0
65 15573.7 217142 0 0 0 0 1 0 0 0 0 0 0
66 17495.0 219452 0 0 0 0 0 1 0 0 0 0 0
67 14155.8 256446 0 0 0 0 0 0 1 0 0 0 0
68 14913.9 265845 0 0 0 0 0 0 0 1 0 0 0
69 17250.4 248624 0 0 0 0 0 0 0 0 1 0 0
70 15879.8 241114 0 0 0 0 0 0 0 0 0 1 0
71 17647.8 229245 0 0 0 0 0 0 0 0 0 0 1
72 17749.9 231805 0 0 0 0 0 0 0 0 0 0 0
73 17111.8 219277 1 0 0 0 0 0 0 0 0 0 0
74 16934.8 219313 0 1 0 0 0 0 0 0 0 0 0
75 20280.0 212610 0 0 1 0 0 0 0 0 0 0 0
76 16238.2 214771 0 0 0 1 0 0 0 0 0 0 0
77 17896.1 211142 0 0 0 0 1 0 0 0 0 0 0
78 18089.3 211457 0 0 0 0 0 1 0 0 0 0 0
79 15660.0 240048 0 0 0 0 0 0 1 0 0 0 0
80 16162.4 240636 0 0 0 0 0 0 0 1 0 0 0
81 17850.1 230580 0 0 0 0 0 0 0 0 1 0 0
82 18520.4 208795 0 0 0 0 0 0 0 0 0 1 0
83 18524.7 197922 0 0 0 0 0 0 0 0 0 0 1
84 16843.7 194596 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) werkloosheid M1 M2 M3
7.372e+03 3.832e-02 -9.334e+02 -4.780e+02 1.564e+03
M4 M5 M6 M7 M8
-2.138e+02 3.223e+02 8.532e+02 -2.374e+03 -3.003e+03
M9 M10 M11
-5.459e+02 2.043e+02 4.135e+02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2306.5 -936.9 -291.1 797.2 3196.6
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.372e+03 1.298e+03 5.681 2.74e-07 ***
werkloosheid 3.832e-02 5.893e-03 6.503 9.46e-09 ***
M1 -9.334e+02 7.380e+02 -1.265 0.210068
M2 -4.780e+02 7.387e+02 -0.647 0.519707
M3 1.564e+03 7.406e+02 2.112 0.038188 *
M4 -2.138e+02 7.423e+02 -0.288 0.774198
M5 3.223e+02 7.452e+02 0.433 0.666679
M6 8.532e+02 7.440e+02 1.147 0.255341
M7 -2.374e+03 7.457e+02 -3.184 0.002155 **
M8 -3.003e+03 7.588e+02 -3.957 0.000178 ***
M9 -5.459e+02 7.433e+02 -0.734 0.465130
M10 2.043e+02 7.389e+02 0.277 0.782949
M11 4.135e+02 7.379e+02 0.560 0.577054
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1380 on 71 degrees of freedom
Multiple R-squared: 0.5118, Adjusted R-squared: 0.4293
F-statistic: 6.202 on 12 and 71 DF, p-value: 2.495e-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,] 7.417361e-02 1.483472e-01 0.9258264
[2,] 3.133534e-02 6.267068e-02 0.9686647
[3,] 1.575005e-02 3.150010e-02 0.9842500
[4,] 5.137455e-03 1.027491e-02 0.9948625
[5,] 6.379659e-03 1.275932e-02 0.9936203
[6,] 2.804851e-03 5.609702e-03 0.9971951
[7,] 2.102031e-03 4.204063e-03 0.9978980
[8,] 8.184457e-04 1.636891e-03 0.9991816
[9,] 2.960303e-04 5.920606e-04 0.9997040
[10,] 1.905264e-04 3.810528e-04 0.9998095
[11,] 7.186624e-05 1.437325e-04 0.9999281
[12,] 2.861291e-05 5.722583e-05 0.9999714
[13,] 1.663205e-04 3.326410e-04 0.9998337
[14,] 6.234632e-05 1.246926e-04 0.9999377
[15,] 3.765031e-05 7.530062e-05 0.9999623
[16,] 1.706452e-04 3.412904e-04 0.9998294
[17,] 7.124934e-05 1.424987e-04 0.9999288
[18,] 5.752893e-05 1.150579e-04 0.9999425
[19,] 8.347985e-05 1.669597e-04 0.9999165
[20,] 4.088050e-05 8.176100e-05 0.9999591
[21,] 4.715997e-05 9.431993e-05 0.9999528
[22,] 5.866022e-05 1.173204e-04 0.9999413
[23,] 3.217203e-05 6.434407e-05 0.9999678
[24,] 2.146253e-05 4.292505e-05 0.9999785
[25,] 1.534868e-05 3.069736e-05 0.9999847
[26,] 1.144299e-05 2.288598e-05 0.9999886
[27,] 1.496236e-05 2.992472e-05 0.9999850
[28,] 1.185622e-05 2.371244e-05 0.9999881
[29,] 1.970258e-05 3.940516e-05 0.9999803
[30,] 3.076725e-05 6.153449e-05 0.9999692
[31,] 2.742697e-05 5.485393e-05 0.9999726
[32,] 1.122500e-04 2.244999e-04 0.9998878
[33,] 2.592676e-04 5.185352e-04 0.9997407
[34,] 4.444766e-04 8.889533e-04 0.9995555
[35,] 7.966530e-04 1.593306e-03 0.9992033
[36,] 4.359701e-03 8.719402e-03 0.9956403
[37,] 7.739142e-03 1.547828e-02 0.9922609
[38,] 6.888190e-02 1.377638e-01 0.9311181
[39,] 2.322400e-01 4.644800e-01 0.7677600
[40,] 2.603460e-01 5.206919e-01 0.7396540
[41,] 3.057920e-01 6.115840e-01 0.6942080
[42,] 3.837075e-01 7.674151e-01 0.6162925
[43,] 3.112103e-01 6.224206e-01 0.6887897
[44,] 3.334323e-01 6.668646e-01 0.6665677
[45,] 2.912310e-01 5.824619e-01 0.7087690
[46,] 3.018203e-01 6.036407e-01 0.6981797
[47,] 2.659989e-01 5.319979e-01 0.7340011
[48,] 4.320239e-01 8.640478e-01 0.5679761
[49,] 3.470824e-01 6.941648e-01 0.6529176
[50,] 5.165545e-01 9.668909e-01 0.4834455
[51,] 4.412690e-01 8.825380e-01 0.5587310
[52,] 3.736257e-01 7.472513e-01 0.6263743
[53,] 2.700494e-01 5.400987e-01 0.7299506
> postscript(file="/var/www/html/rcomp/tmp/189su1229262244.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/2v9yy1229262244.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/3j1je1229262244.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/44fte1229262244.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/5qurt1229262244.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
-1787.10059 -366.70691 -485.16025 -1198.78247 800.54818 -862.11478
7 8 9 10 11 12
38.37236 590.71429 445.45375 610.37876 1536.32545 796.01887
13 14 15 16 17 18
1797.03755 1207.82567 663.83461 717.78092 1170.02719 384.83897
19 20 21 22 23 24
629.73158 292.61558 -981.34194 -723.32545 -1047.69370 -1811.27703
25 26 27 28 29 30
-426.76851 -652.93917 -1293.57843 184.80122 -302.48973 -1049.06411
31 32 33 34 35 36
-102.12402 -897.74276 -720.61569 17.39645 -1532.09135 -832.81403
37 38 39 40 41 42
-93.37126 -565.57655 -1485.92016 -570.62243 -1634.84256 -1871.49446
43 44 45 46 47 48
-939.62084 -2306.46081 -1548.97233 -929.93784 -2197.27579 -1283.51242
49 50 51 52 53 54
-1171.21270 -1071.28922 -427.04742 -255.77544 -1702.71017 776.27778
55 56 57 58 59 60
-421.47113 -608.60660 -279.76227 -981.62008 -991.92038 -378.43253
61 62 63 64 65 66
-589.16039 -188.05910 -168.76076 272.60684 -441.43005 860.43931
67 68 69 70 71 72
-668.73949 357.47059 897.03537 -935.99735 1077.71041 1495.16180
73 74 75 76 77 78
2270.57589 1636.74527 3196.63241 849.99137 2110.89715 1761.11731
79 80 81 82 83 84
1463.85154 2572.00972 2188.20311 2943.10551 3154.94536 2014.85533
> postscript(file="/var/www/html/rcomp/tmp/6ef5n1229262244.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 -1787.10059 NA
1 -366.70691 -1787.10059
2 -485.16025 -366.70691
3 -1198.78247 -485.16025
4 800.54818 -1198.78247
5 -862.11478 800.54818
6 38.37236 -862.11478
7 590.71429 38.37236
8 445.45375 590.71429
9 610.37876 445.45375
10 1536.32545 610.37876
11 796.01887 1536.32545
12 1797.03755 796.01887
13 1207.82567 1797.03755
14 663.83461 1207.82567
15 717.78092 663.83461
16 1170.02719 717.78092
17 384.83897 1170.02719
18 629.73158 384.83897
19 292.61558 629.73158
20 -981.34194 292.61558
21 -723.32545 -981.34194
22 -1047.69370 -723.32545
23 -1811.27703 -1047.69370
24 -426.76851 -1811.27703
25 -652.93917 -426.76851
26 -1293.57843 -652.93917
27 184.80122 -1293.57843
28 -302.48973 184.80122
29 -1049.06411 -302.48973
30 -102.12402 -1049.06411
31 -897.74276 -102.12402
32 -720.61569 -897.74276
33 17.39645 -720.61569
34 -1532.09135 17.39645
35 -832.81403 -1532.09135
36 -93.37126 -832.81403
37 -565.57655 -93.37126
38 -1485.92016 -565.57655
39 -570.62243 -1485.92016
40 -1634.84256 -570.62243
41 -1871.49446 -1634.84256
42 -939.62084 -1871.49446
43 -2306.46081 -939.62084
44 -1548.97233 -2306.46081
45 -929.93784 -1548.97233
46 -2197.27579 -929.93784
47 -1283.51242 -2197.27579
48 -1171.21270 -1283.51242
49 -1071.28922 -1171.21270
50 -427.04742 -1071.28922
51 -255.77544 -427.04742
52 -1702.71017 -255.77544
53 776.27778 -1702.71017
54 -421.47113 776.27778
55 -608.60660 -421.47113
56 -279.76227 -608.60660
57 -981.62008 -279.76227
58 -991.92038 -981.62008
59 -378.43253 -991.92038
60 -589.16039 -378.43253
61 -188.05910 -589.16039
62 -168.76076 -188.05910
63 272.60684 -168.76076
64 -441.43005 272.60684
65 860.43931 -441.43005
66 -668.73949 860.43931
67 357.47059 -668.73949
68 897.03537 357.47059
69 -935.99735 897.03537
70 1077.71041 -935.99735
71 1495.16180 1077.71041
72 2270.57589 1495.16180
73 1636.74527 2270.57589
74 3196.63241 1636.74527
75 849.99137 3196.63241
76 2110.89715 849.99137
77 1761.11731 2110.89715
78 1463.85154 1761.11731
79 2572.00972 1463.85154
80 2188.20311 2572.00972
81 2943.10551 2188.20311
82 3154.94536 2943.10551
83 2014.85533 3154.94536
84 NA 2014.85533
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -366.70691 -1787.10059
[2,] -485.16025 -366.70691
[3,] -1198.78247 -485.16025
[4,] 800.54818 -1198.78247
[5,] -862.11478 800.54818
[6,] 38.37236 -862.11478
[7,] 590.71429 38.37236
[8,] 445.45375 590.71429
[9,] 610.37876 445.45375
[10,] 1536.32545 610.37876
[11,] 796.01887 1536.32545
[12,] 1797.03755 796.01887
[13,] 1207.82567 1797.03755
[14,] 663.83461 1207.82567
[15,] 717.78092 663.83461
[16,] 1170.02719 717.78092
[17,] 384.83897 1170.02719
[18,] 629.73158 384.83897
[19,] 292.61558 629.73158
[20,] -981.34194 292.61558
[21,] -723.32545 -981.34194
[22,] -1047.69370 -723.32545
[23,] -1811.27703 -1047.69370
[24,] -426.76851 -1811.27703
[25,] -652.93917 -426.76851
[26,] -1293.57843 -652.93917
[27,] 184.80122 -1293.57843
[28,] -302.48973 184.80122
[29,] -1049.06411 -302.48973
[30,] -102.12402 -1049.06411
[31,] -897.74276 -102.12402
[32,] -720.61569 -897.74276
[33,] 17.39645 -720.61569
[34,] -1532.09135 17.39645
[35,] -832.81403 -1532.09135
[36,] -93.37126 -832.81403
[37,] -565.57655 -93.37126
[38,] -1485.92016 -565.57655
[39,] -570.62243 -1485.92016
[40,] -1634.84256 -570.62243
[41,] -1871.49446 -1634.84256
[42,] -939.62084 -1871.49446
[43,] -2306.46081 -939.62084
[44,] -1548.97233 -2306.46081
[45,] -929.93784 -1548.97233
[46,] -2197.27579 -929.93784
[47,] -1283.51242 -2197.27579
[48,] -1171.21270 -1283.51242
[49,] -1071.28922 -1171.21270
[50,] -427.04742 -1071.28922
[51,] -255.77544 -427.04742
[52,] -1702.71017 -255.77544
[53,] 776.27778 -1702.71017
[54,] -421.47113 776.27778
[55,] -608.60660 -421.47113
[56,] -279.76227 -608.60660
[57,] -981.62008 -279.76227
[58,] -991.92038 -981.62008
[59,] -378.43253 -991.92038
[60,] -589.16039 -378.43253
[61,] -188.05910 -589.16039
[62,] -168.76076 -188.05910
[63,] 272.60684 -168.76076
[64,] -441.43005 272.60684
[65,] 860.43931 -441.43005
[66,] -668.73949 860.43931
[67,] 357.47059 -668.73949
[68,] 897.03537 357.47059
[69,] -935.99735 897.03537
[70,] 1077.71041 -935.99735
[71,] 1495.16180 1077.71041
[72,] 2270.57589 1495.16180
[73,] 1636.74527 2270.57589
[74,] 3196.63241 1636.74527
[75,] 849.99137 3196.63241
[76,] 2110.89715 849.99137
[77,] 1761.11731 2110.89715
[78,] 1463.85154 1761.11731
[79,] 2572.00972 1463.85154
[80,] 2188.20311 2572.00972
[81,] 2943.10551 2188.20311
[82,] 3154.94536 2943.10551
[83,] 2014.85533 3154.94536
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -366.70691 -1787.10059
2 -485.16025 -366.70691
3 -1198.78247 -485.16025
4 800.54818 -1198.78247
5 -862.11478 800.54818
6 38.37236 -862.11478
7 590.71429 38.37236
8 445.45375 590.71429
9 610.37876 445.45375
10 1536.32545 610.37876
11 796.01887 1536.32545
12 1797.03755 796.01887
13 1207.82567 1797.03755
14 663.83461 1207.82567
15 717.78092 663.83461
16 1170.02719 717.78092
17 384.83897 1170.02719
18 629.73158 384.83897
19 292.61558 629.73158
20 -981.34194 292.61558
21 -723.32545 -981.34194
22 -1047.69370 -723.32545
23 -1811.27703 -1047.69370
24 -426.76851 -1811.27703
25 -652.93917 -426.76851
26 -1293.57843 -652.93917
27 184.80122 -1293.57843
28 -302.48973 184.80122
29 -1049.06411 -302.48973
30 -102.12402 -1049.06411
31 -897.74276 -102.12402
32 -720.61569 -897.74276
33 17.39645 -720.61569
34 -1532.09135 17.39645
35 -832.81403 -1532.09135
36 -93.37126 -832.81403
37 -565.57655 -93.37126
38 -1485.92016 -565.57655
39 -570.62243 -1485.92016
40 -1634.84256 -570.62243
41 -1871.49446 -1634.84256
42 -939.62084 -1871.49446
43 -2306.46081 -939.62084
44 -1548.97233 -2306.46081
45 -929.93784 -1548.97233
46 -2197.27579 -929.93784
47 -1283.51242 -2197.27579
48 -1171.21270 -1283.51242
49 -1071.28922 -1171.21270
50 -427.04742 -1071.28922
51 -255.77544 -427.04742
52 -1702.71017 -255.77544
53 776.27778 -1702.71017
54 -421.47113 776.27778
55 -608.60660 -421.47113
56 -279.76227 -608.60660
57 -981.62008 -279.76227
58 -991.92038 -981.62008
59 -378.43253 -991.92038
60 -589.16039 -378.43253
61 -188.05910 -589.16039
62 -168.76076 -188.05910
63 272.60684 -168.76076
64 -441.43005 272.60684
65 860.43931 -441.43005
66 -668.73949 860.43931
67 357.47059 -668.73949
68 897.03537 357.47059
69 -935.99735 897.03537
70 1077.71041 -935.99735
71 1495.16180 1077.71041
72 2270.57589 1495.16180
73 1636.74527 2270.57589
74 3196.63241 1636.74527
75 849.99137 3196.63241
76 2110.89715 849.99137
77 1761.11731 2110.89715
78 1463.85154 1761.11731
79 2572.00972 1463.85154
80 2188.20311 2572.00972
81 2943.10551 2188.20311
82 3154.94536 2943.10551
83 2014.85533 3154.94536
> 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/7f00f1229262244.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/8a1vm1229262244.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/9fjfq1229262244.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/10z3311229262244.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/11mjzo1229262244.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/12mizr1229262244.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/13lius1229262244.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/14v74d1229262244.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/15qtte1229262244.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/16zno51229262244.tab")
+ }
>
> system("convert tmp/189su1229262244.ps tmp/189su1229262244.png")
> system("convert tmp/2v9yy1229262244.ps tmp/2v9yy1229262244.png")
> system("convert tmp/3j1je1229262244.ps tmp/3j1je1229262244.png")
> system("convert tmp/44fte1229262244.ps tmp/44fte1229262244.png")
> system("convert tmp/5qurt1229262244.ps tmp/5qurt1229262244.png")
> system("convert tmp/6ef5n1229262244.ps tmp/6ef5n1229262244.png")
> system("convert tmp/7f00f1229262244.ps tmp/7f00f1229262244.png")
> system("convert tmp/8a1vm1229262244.ps tmp/8a1vm1229262244.png")
> system("convert tmp/9fjfq1229262244.ps tmp/9fjfq1229262244.png")
> system("convert tmp/10z3311229262244.ps tmp/10z3311229262244.png")
>
>
> proc.time()
user system elapsed
2.641 1.577 3.149