R version 2.7.0 (2008-04-22)
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
+ ,7.5
+ ,13182.1
+ ,7.2
+ ,14800.1
+ ,6.9
+ ,12150.7
+ ,6.7
+ ,14478.2
+ ,6.4
+ ,13253.9
+ ,6.3
+ ,12036.8
+ ,6.8
+ ,12653.2
+ ,7.3
+ ,14035.4
+ ,7.1
+ ,14571.4
+ ,7.1
+ ,15400.9
+ ,6.8
+ ,14283.2
+ ,6.5
+ ,14485.3
+ ,6.3
+ ,14196.3
+ ,6.1
+ ,15559.1
+ ,6.1
+ ,13767.4
+ ,6.3
+ ,14634
+ ,6.3
+ ,14381.1
+ ,6
+ ,12509.9
+ ,6.2
+ ,12122.3
+ ,6.4
+ ,13122.3
+ ,6.8
+ ,13908.7
+ ,7.5
+ ,13456.5
+ ,7.5
+ ,12441.6
+ ,7.6
+ ,12953
+ ,7.6
+ ,13057.2
+ ,7.4
+ ,14350.1
+ ,7.3
+ ,13830.2
+ ,7.1
+ ,13755.5
+ ,6.9
+ ,13574.4
+ ,6.8
+ ,12802.6
+ ,7.5
+ ,11737.3
+ ,7.6
+ ,13850.2
+ ,7.8
+ ,15081.8
+ ,8
+ ,13653.3
+ ,8.1
+ ,14019.1
+ ,8.2
+ ,13962
+ ,8.3
+ ,13768.7
+ ,8.2
+ ,14747.1
+ ,8
+ ,13858.1
+ ,7.9
+ ,13188
+ ,7.6
+ ,13693.1
+ ,7.6
+ ,12970
+ ,8.2
+ ,11392.8
+ ,8.3
+ ,13985.2
+ ,8.4
+ ,14994.7
+ ,8.4
+ ,13584.7
+ ,8.4
+ ,14257.8
+ ,8.6
+ ,13553.4
+ ,8.9
+ ,14007.3
+ ,8.8
+ ,16535.8
+ ,8.3
+ ,14721.4
+ ,7.5
+ ,13664.6
+ ,7.2
+ ,16805.9
+ ,7.5
+ ,13829.4
+ ,8.8
+ ,13735.6
+ ,9.3
+ ,15870.5
+ ,9.3
+ ,15962.4
+ ,8.7
+ ,15744.1
+ ,8.2
+ ,16083.7
+ ,8.3
+ ,14863.9
+ ,8.5
+ ,15533.1
+ ,8.6
+ ,17473.1
+ ,8.6
+ ,15925.5
+ ,8.2
+ ,15573.7
+ ,8.1
+ ,17495
+ ,8
+ ,14155.8
+ ,8.6
+ ,14913.9
+ ,8.7
+ ,17250.4
+ ,8.8
+ ,15879.8
+ ,8.5
+ ,17647.8
+ ,8.4
+ ,17749.9
+ ,8.5
+ ,17111.8
+ ,8.7
+ ,16934.8
+ ,8.7
+ ,20280
+ ,8.6
+ ,16238.2
+ ,8.5
+ ,17896.1
+ ,8.3
+ ,18089.3
+ ,8.1
+ ,15660
+ ,8.2
+ ,16162.4
+ ,8.1
+ ,17850.1
+ ,8.1
+ ,18520.4
+ ,7.9
+ ,18524.7
+ ,7.9
+ ,16843.7
+ ,7.9)
+ ,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 7.5 1 0 0 0 0 0 0 0 0 0 0
2 13182.1 7.2 0 1 0 0 0 0 0 0 0 0 0
3 14800.1 6.9 0 0 1 0 0 0 0 0 0 0 0
4 12150.7 6.7 0 0 0 1 0 0 0 0 0 0 0
5 14478.2 6.4 0 0 0 0 1 0 0 0 0 0 0
6 13253.9 6.3 0 0 0 0 0 1 0 0 0 0 0
7 12036.8 6.8 0 0 0 0 0 0 1 0 0 0 0
8 12653.2 7.3 0 0 0 0 0 0 0 1 0 0 0
9 14035.4 7.1 0 0 0 0 0 0 0 0 1 0 0
10 14571.4 7.1 0 0 0 0 0 0 0 0 0 1 0
11 15400.9 6.8 0 0 0 0 0 0 0 0 0 0 1
12 14283.2 6.5 0 0 0 0 0 0 0 0 0 0 0
13 14485.3 6.3 1 0 0 0 0 0 0 0 0 0 0
14 14196.3 6.1 0 1 0 0 0 0 0 0 0 0 0
15 15559.1 6.1 0 0 1 0 0 0 0 0 0 0 0
16 13767.4 6.3 0 0 0 1 0 0 0 0 0 0 0
17 14634.0 6.3 0 0 0 0 1 0 0 0 0 0 0
18 14381.1 6.0 0 0 0 0 0 1 0 0 0 0 0
19 12509.9 6.2 0 0 0 0 0 0 1 0 0 0 0
20 12122.3 6.4 0 0 0 0 0 0 0 1 0 0 0
21 13122.3 6.8 0 0 0 0 0 0 0 0 1 0 0
22 13908.7 7.5 0 0 0 0 0 0 0 0 0 1 0
23 13456.5 7.5 0 0 0 0 0 0 0 0 0 0 1
24 12441.6 7.6 0 0 0 0 0 0 0 0 0 0 0
25 12953.0 7.6 1 0 0 0 0 0 0 0 0 0 0
26 13057.2 7.4 0 1 0 0 0 0 0 0 0 0 0
27 14350.1 7.3 0 0 1 0 0 0 0 0 0 0 0
28 13830.2 7.1 0 0 0 1 0 0 0 0 0 0 0
29 13755.5 6.9 0 0 0 0 1 0 0 0 0 0 0
30 13574.4 6.8 0 0 0 0 0 1 0 0 0 0 0
31 12802.6 7.5 0 0 0 0 0 0 1 0 0 0 0
32 11737.3 7.6 0 0 0 0 0 0 0 1 0 0 0
33 13850.2 7.8 0 0 0 0 0 0 0 0 1 0 0
34 15081.8 8.0 0 0 0 0 0 0 0 0 0 1 0
35 13653.3 8.1 0 0 0 0 0 0 0 0 0 0 1
36 14019.1 8.2 0 0 0 0 0 0 0 0 0 0 0
37 13962.0 8.3 1 0 0 0 0 0 0 0 0 0 0
38 13768.7 8.2 0 1 0 0 0 0 0 0 0 0 0
39 14747.1 8.0 0 0 1 0 0 0 0 0 0 0 0
40 13858.1 7.9 0 0 0 1 0 0 0 0 0 0 0
41 13188.0 7.6 0 0 0 0 1 0 0 0 0 0 0
42 13693.1 7.6 0 0 0 0 0 1 0 0 0 0 0
43 12970.0 8.2 0 0 0 0 0 0 1 0 0 0 0
44 11392.8 8.3 0 0 0 0 0 0 0 1 0 0 0
45 13985.2 8.4 0 0 0 0 0 0 0 0 1 0 0
46 14994.7 8.4 0 0 0 0 0 0 0 0 0 1 0
47 13584.7 8.4 0 0 0 0 0 0 0 0 0 0 1
48 14257.8 8.6 0 0 0 0 0 0 0 0 0 0 0
49 13553.4 8.9 1 0 0 0 0 0 0 0 0 0 0
50 14007.3 8.8 0 1 0 0 0 0 0 0 0 0 0
51 16535.8 8.3 0 0 1 0 0 0 0 0 0 0 0
52 14721.4 7.5 0 0 0 1 0 0 0 0 0 0 0
53 13664.6 7.2 0 0 0 0 1 0 0 0 0 0 0
54 16805.9 7.5 0 0 0 0 0 1 0 0 0 0 0
55 13829.4 8.8 0 0 0 0 0 0 1 0 0 0 0
56 13735.6 9.3 0 0 0 0 0 0 0 1 0 0 0
57 15870.5 9.3 0 0 0 0 0 0 0 0 1 0 0
58 15962.4 8.7 0 0 0 0 0 0 0 0 0 1 0
59 15744.1 8.2 0 0 0 0 0 0 0 0 0 0 1
60 16083.7 8.3 0 0 0 0 0 0 0 0 0 0 0
61 14863.9 8.5 1 0 0 0 0 0 0 0 0 0 0
62 15533.1 8.6 0 1 0 0 0 0 0 0 0 0 0
63 17473.1 8.6 0 0 1 0 0 0 0 0 0 0 0
64 15925.5 8.2 0 0 0 1 0 0 0 0 0 0 0
65 15573.7 8.1 0 0 0 0 1 0 0 0 0 0 0
66 17495.0 8.0 0 0 0 0 0 1 0 0 0 0 0
67 14155.8 8.6 0 0 0 0 0 0 1 0 0 0 0
68 14913.9 8.7 0 0 0 0 0 0 0 1 0 0 0
69 17250.4 8.8 0 0 0 0 0 0 0 0 1 0 0
70 15879.8 8.5 0 0 0 0 0 0 0 0 0 1 0
71 17647.8 8.4 0 0 0 0 0 0 0 0 0 0 1
72 17749.9 8.5 0 0 0 0 0 0 0 0 0 0 0
73 17111.8 8.7 1 0 0 0 0 0 0 0 0 0 0
74 16934.8 8.7 0 1 0 0 0 0 0 0 0 0 0
75 20280.0 8.6 0 0 1 0 0 0 0 0 0 0 0
76 16238.2 8.5 0 0 0 1 0 0 0 0 0 0 0
77 17896.1 8.3 0 0 0 0 1 0 0 0 0 0 0
78 18089.3 8.1 0 0 0 0 0 1 0 0 0 0 0
79 15660.0 8.2 0 0 0 0 0 0 1 0 0 0 0
80 16162.4 8.1 0 0 0 0 0 0 0 1 0 0 0
81 17850.1 8.1 0 0 0 0 0 0 0 0 1 0 0
82 18520.4 7.9 0 0 0 0 0 0 0 0 0 1 0
83 18524.7 7.9 0 0 0 0 0 0 0 0 0 0 1
84 16843.7 7.9 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
6897.17 1032.35 -1057.37 -625.73 1417.79
M4 M5 M6 M7 M8
-239.64 352.34 1012.17 -1481.78 -1866.39
M9 M10 M11
-62.51 389.15 377.53
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2361.8 -1022.6 -169.3 885.9 3094.4
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6897.17 1771.68 3.893 0.000221 ***
Werkloosheid 1032.35 211.19 4.888 6.1e-06 ***
M1 -1057.37 806.19 -1.312 0.193896
M2 -625.73 806.38 -0.776 0.440341
M3 1417.79 808.00 1.755 0.083626 .
M4 -239.64 812.67 -0.295 0.768944
M5 352.34 819.08 0.430 0.668375
M6 1012.17 821.88 1.232 0.222188
M7 -1481.78 807.13 -1.836 0.070564 .
M8 -1866.39 806.18 -2.315 0.023502 *
M9 -62.51 806.45 -0.078 0.938437
M10 389.15 806.31 0.483 0.630847
M11 377.53 806.22 0.468 0.641026
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1508 on 71 degrees of freedom
Multiple R-squared: 0.4172, Adjusted R-squared: 0.3186
F-statistic: 4.235 on 12 and 71 DF, p-value: 5.162e-05
> 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,] 9.167459e-02 1.833492e-01 0.9083254
[2,] 3.165864e-02 6.331728e-02 0.9683414
[3,] 1.294788e-02 2.589577e-02 0.9870521
[4,] 5.223812e-03 1.044762e-02 0.9947762
[5,] 1.372358e-02 2.744717e-02 0.9862764
[6,] 1.085748e-02 2.171496e-02 0.9891425
[7,] 4.456813e-03 8.913626e-03 0.9955432
[8,] 2.896023e-03 5.792046e-03 0.9971040
[9,] 1.342332e-03 2.684664e-03 0.9986577
[10,] 8.296275e-04 1.659255e-03 0.9991704
[11,] 3.305939e-04 6.611878e-04 0.9996694
[12,] 1.284516e-04 2.569031e-04 0.9998715
[13,] 2.300372e-04 4.600744e-04 0.9997700
[14,] 8.828796e-05 1.765759e-04 0.9999117
[15,] 4.136229e-05 8.272459e-05 0.9999586
[16,] 5.921537e-05 1.184307e-04 0.9999408
[17,] 2.365550e-05 4.731100e-05 0.9999763
[18,] 1.743853e-05 3.487707e-05 0.9999826
[19,] 1.904585e-05 3.809169e-05 0.9999810
[20,] 1.085908e-05 2.171816e-05 0.9999891
[21,] 1.529027e-05 3.058055e-05 0.9999847
[22,] 2.255664e-05 4.511328e-05 0.9999774
[23,] 1.413948e-05 2.827896e-05 0.9999859
[24,] 1.159199e-05 2.318397e-05 0.9999884
[25,] 9.262500e-06 1.852500e-05 0.9999907
[26,] 8.265092e-06 1.653018e-05 0.9999917
[27,] 1.225652e-05 2.451304e-05 0.9999877
[28,] 9.463241e-06 1.892648e-05 0.9999905
[29,] 1.788579e-05 3.577159e-05 0.9999821
[30,] 2.991751e-05 5.983501e-05 0.9999701
[31,] 2.409239e-05 4.818479e-05 0.9999759
[32,] 7.182963e-05 1.436593e-04 0.9999282
[33,] 1.293650e-04 2.587301e-04 0.9998706
[34,] 1.663322e-04 3.326645e-04 0.9998337
[35,] 2.155106e-04 4.310213e-04 0.9997845
[36,] 9.422465e-04 1.884493e-03 0.9990578
[37,] 1.407898e-03 2.815795e-03 0.9985921
[38,] 1.466439e-02 2.932879e-02 0.9853356
[39,] 8.748402e-02 1.749680e-01 0.9125160
[40,] 7.244284e-02 1.448857e-01 0.9275572
[41,] 6.603541e-02 1.320708e-01 0.9339646
[42,] 6.497696e-02 1.299539e-01 0.9350230
[43,] 5.164211e-02 1.032842e-01 0.9483579
[44,] 1.042528e-01 2.085056e-01 0.8957472
[45,] 1.252475e-01 2.504951e-01 0.8747525
[46,] 2.006766e-01 4.013533e-01 0.7993234
[47,] 2.131065e-01 4.262131e-01 0.7868935
[48,] 4.435407e-01 8.870814e-01 0.5564593
[49,] 3.918540e-01 7.837081e-01 0.6081460
[50,] 6.684364e-01 6.631273e-01 0.3315636
[51,] 6.288134e-01 7.423732e-01 0.3711866
[52,] 5.385864e-01 9.228273e-01 0.4614136
[53,] 4.085335e-01 8.170671e-01 0.5914665
> postscript(file="/var/www/html/rcomp/tmp/12z6o1228659710.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/21x2c1228659710.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/3vt9i1228659710.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/4brwd1228659710.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/5n8fu1228659710.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
-2027.948097 -522.282629 -638.094448 -1423.590234 621.630732 -1159.259196
7 8 9 10 11 12
-398.591158 86.245942 -128.967623 -44.620552 1106.202059 675.737209
13 14 15 16 17 18
2141.674928 1627.505144 946.787569 606.050774 880.665984 277.646561
19 20 21 22 23 24
693.920354 484.463211 -732.361867 -1120.261561 -1560.844706 -2301.450564
25 26 27 28 29 30
-732.683349 -853.653133 -1501.035456 -157.031242 -617.245528 -1354.935456
31 32 33 34 35 36
-355.437923 -1139.359814 -1036.814388 -463.337821 -1983.456218 -1343.362077
37 38 39 40 41 42
-446.330114 -968.035150 -1826.682221 -955.013259 -1907.392293 -2062.117473
43 44 45 46 47 48
-910.684688 -2206.506579 -1521.225900 -963.378830 -2361.761975 -1517.603085
49 50 51 52 53 54
-1474.341627 -1348.846663 -347.687977 321.227749 -1017.851285 1153.917779
55 56 57 58 59 60
-670.696200 -896.059100 -565.043169 -305.384586 4.108529 618.002671
61 62 63 64 65 66
249.099382 383.423842 279.906267 802.680984 -37.868553 1326.841519
67 68 69 70 71 72
-137.825696 901.652413 1331.033091 -181.514082 1701.338025 2077.732167
73 74 75 76 77 78
2290.528878 1681.888589 3086.806267 805.675228 2078.060942 1817.906267
79 80 81 82 83 84
1779.315312 2769.563926 2653.379856 3078.497431 3094.414286 1790.943679
> postscript(file="/var/www/html/rcomp/tmp/6udyv1228659710.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 -2027.948097 NA
1 -522.282629 -2027.948097
2 -638.094448 -522.282629
3 -1423.590234 -638.094448
4 621.630732 -1423.590234
5 -1159.259196 621.630732
6 -398.591158 -1159.259196
7 86.245942 -398.591158
8 -128.967623 86.245942
9 -44.620552 -128.967623
10 1106.202059 -44.620552
11 675.737209 1106.202059
12 2141.674928 675.737209
13 1627.505144 2141.674928
14 946.787569 1627.505144
15 606.050774 946.787569
16 880.665984 606.050774
17 277.646561 880.665984
18 693.920354 277.646561
19 484.463211 693.920354
20 -732.361867 484.463211
21 -1120.261561 -732.361867
22 -1560.844706 -1120.261561
23 -2301.450564 -1560.844706
24 -732.683349 -2301.450564
25 -853.653133 -732.683349
26 -1501.035456 -853.653133
27 -157.031242 -1501.035456
28 -617.245528 -157.031242
29 -1354.935456 -617.245528
30 -355.437923 -1354.935456
31 -1139.359814 -355.437923
32 -1036.814388 -1139.359814
33 -463.337821 -1036.814388
34 -1983.456218 -463.337821
35 -1343.362077 -1983.456218
36 -446.330114 -1343.362077
37 -968.035150 -446.330114
38 -1826.682221 -968.035150
39 -955.013259 -1826.682221
40 -1907.392293 -955.013259
41 -2062.117473 -1907.392293
42 -910.684688 -2062.117473
43 -2206.506579 -910.684688
44 -1521.225900 -2206.506579
45 -963.378830 -1521.225900
46 -2361.761975 -963.378830
47 -1517.603085 -2361.761975
48 -1474.341627 -1517.603085
49 -1348.846663 -1474.341627
50 -347.687977 -1348.846663
51 321.227749 -347.687977
52 -1017.851285 321.227749
53 1153.917779 -1017.851285
54 -670.696200 1153.917779
55 -896.059100 -670.696200
56 -565.043169 -896.059100
57 -305.384586 -565.043169
58 4.108529 -305.384586
59 618.002671 4.108529
60 249.099382 618.002671
61 383.423842 249.099382
62 279.906267 383.423842
63 802.680984 279.906267
64 -37.868553 802.680984
65 1326.841519 -37.868553
66 -137.825696 1326.841519
67 901.652413 -137.825696
68 1331.033091 901.652413
69 -181.514082 1331.033091
70 1701.338025 -181.514082
71 2077.732167 1701.338025
72 2290.528878 2077.732167
73 1681.888589 2290.528878
74 3086.806267 1681.888589
75 805.675228 3086.806267
76 2078.060942 805.675228
77 1817.906267 2078.060942
78 1779.315312 1817.906267
79 2769.563926 1779.315312
80 2653.379856 2769.563926
81 3078.497431 2653.379856
82 3094.414286 3078.497431
83 1790.943679 3094.414286
84 NA 1790.943679
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -522.282629 -2027.948097
[2,] -638.094448 -522.282629
[3,] -1423.590234 -638.094448
[4,] 621.630732 -1423.590234
[5,] -1159.259196 621.630732
[6,] -398.591158 -1159.259196
[7,] 86.245942 -398.591158
[8,] -128.967623 86.245942
[9,] -44.620552 -128.967623
[10,] 1106.202059 -44.620552
[11,] 675.737209 1106.202059
[12,] 2141.674928 675.737209
[13,] 1627.505144 2141.674928
[14,] 946.787569 1627.505144
[15,] 606.050774 946.787569
[16,] 880.665984 606.050774
[17,] 277.646561 880.665984
[18,] 693.920354 277.646561
[19,] 484.463211 693.920354
[20,] -732.361867 484.463211
[21,] -1120.261561 -732.361867
[22,] -1560.844706 -1120.261561
[23,] -2301.450564 -1560.844706
[24,] -732.683349 -2301.450564
[25,] -853.653133 -732.683349
[26,] -1501.035456 -853.653133
[27,] -157.031242 -1501.035456
[28,] -617.245528 -157.031242
[29,] -1354.935456 -617.245528
[30,] -355.437923 -1354.935456
[31,] -1139.359814 -355.437923
[32,] -1036.814388 -1139.359814
[33,] -463.337821 -1036.814388
[34,] -1983.456218 -463.337821
[35,] -1343.362077 -1983.456218
[36,] -446.330114 -1343.362077
[37,] -968.035150 -446.330114
[38,] -1826.682221 -968.035150
[39,] -955.013259 -1826.682221
[40,] -1907.392293 -955.013259
[41,] -2062.117473 -1907.392293
[42,] -910.684688 -2062.117473
[43,] -2206.506579 -910.684688
[44,] -1521.225900 -2206.506579
[45,] -963.378830 -1521.225900
[46,] -2361.761975 -963.378830
[47,] -1517.603085 -2361.761975
[48,] -1474.341627 -1517.603085
[49,] -1348.846663 -1474.341627
[50,] -347.687977 -1348.846663
[51,] 321.227749 -347.687977
[52,] -1017.851285 321.227749
[53,] 1153.917779 -1017.851285
[54,] -670.696200 1153.917779
[55,] -896.059100 -670.696200
[56,] -565.043169 -896.059100
[57,] -305.384586 -565.043169
[58,] 4.108529 -305.384586
[59,] 618.002671 4.108529
[60,] 249.099382 618.002671
[61,] 383.423842 249.099382
[62,] 279.906267 383.423842
[63,] 802.680984 279.906267
[64,] -37.868553 802.680984
[65,] 1326.841519 -37.868553
[66,] -137.825696 1326.841519
[67,] 901.652413 -137.825696
[68,] 1331.033091 901.652413
[69,] -181.514082 1331.033091
[70,] 1701.338025 -181.514082
[71,] 2077.732167 1701.338025
[72,] 2290.528878 2077.732167
[73,] 1681.888589 2290.528878
[74,] 3086.806267 1681.888589
[75,] 805.675228 3086.806267
[76,] 2078.060942 805.675228
[77,] 1817.906267 2078.060942
[78,] 1779.315312 1817.906267
[79,] 2769.563926 1779.315312
[80,] 2653.379856 2769.563926
[81,] 3078.497431 2653.379856
[82,] 3094.414286 3078.497431
[83,] 1790.943679 3094.414286
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -522.282629 -2027.948097
2 -638.094448 -522.282629
3 -1423.590234 -638.094448
4 621.630732 -1423.590234
5 -1159.259196 621.630732
6 -398.591158 -1159.259196
7 86.245942 -398.591158
8 -128.967623 86.245942
9 -44.620552 -128.967623
10 1106.202059 -44.620552
11 675.737209 1106.202059
12 2141.674928 675.737209
13 1627.505144 2141.674928
14 946.787569 1627.505144
15 606.050774 946.787569
16 880.665984 606.050774
17 277.646561 880.665984
18 693.920354 277.646561
19 484.463211 693.920354
20 -732.361867 484.463211
21 -1120.261561 -732.361867
22 -1560.844706 -1120.261561
23 -2301.450564 -1560.844706
24 -732.683349 -2301.450564
25 -853.653133 -732.683349
26 -1501.035456 -853.653133
27 -157.031242 -1501.035456
28 -617.245528 -157.031242
29 -1354.935456 -617.245528
30 -355.437923 -1354.935456
31 -1139.359814 -355.437923
32 -1036.814388 -1139.359814
33 -463.337821 -1036.814388
34 -1983.456218 -463.337821
35 -1343.362077 -1983.456218
36 -446.330114 -1343.362077
37 -968.035150 -446.330114
38 -1826.682221 -968.035150
39 -955.013259 -1826.682221
40 -1907.392293 -955.013259
41 -2062.117473 -1907.392293
42 -910.684688 -2062.117473
43 -2206.506579 -910.684688
44 -1521.225900 -2206.506579
45 -963.378830 -1521.225900
46 -2361.761975 -963.378830
47 -1517.603085 -2361.761975
48 -1474.341627 -1517.603085
49 -1348.846663 -1474.341627
50 -347.687977 -1348.846663
51 321.227749 -347.687977
52 -1017.851285 321.227749
53 1153.917779 -1017.851285
54 -670.696200 1153.917779
55 -896.059100 -670.696200
56 -565.043169 -896.059100
57 -305.384586 -565.043169
58 4.108529 -305.384586
59 618.002671 4.108529
60 249.099382 618.002671
61 383.423842 249.099382
62 279.906267 383.423842
63 802.680984 279.906267
64 -37.868553 802.680984
65 1326.841519 -37.868553
66 -137.825696 1326.841519
67 901.652413 -137.825696
68 1331.033091 901.652413
69 -181.514082 1331.033091
70 1701.338025 -181.514082
71 2077.732167 1701.338025
72 2290.528878 2077.732167
73 1681.888589 2290.528878
74 3086.806267 1681.888589
75 805.675228 3086.806267
76 2078.060942 805.675228
77 1817.906267 2078.060942
78 1779.315312 1817.906267
79 2769.563926 1779.315312
80 2653.379856 2769.563926
81 3078.497431 2653.379856
82 3094.414286 3078.497431
83 1790.943679 3094.414286
> 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/72ffs1228659710.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/89xt91228659710.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/9lc0g1228659710.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/105lhm1228659710.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/113uie1228659710.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/129g711228659710.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/13fvud1228659710.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/14apzv1228659710.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/15aybg1228659710.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/162hju1228659710.tab")
+ }
>
> system("convert tmp/12z6o1228659710.ps tmp/12z6o1228659710.png")
> system("convert tmp/21x2c1228659710.ps tmp/21x2c1228659710.png")
> system("convert tmp/3vt9i1228659710.ps tmp/3vt9i1228659710.png")
> system("convert tmp/4brwd1228659710.ps tmp/4brwd1228659710.png")
> system("convert tmp/5n8fu1228659710.ps tmp/5n8fu1228659710.png")
> system("convert tmp/6udyv1228659710.ps tmp/6udyv1228659710.png")
> system("convert tmp/72ffs1228659710.ps tmp/72ffs1228659710.png")
> system("convert tmp/89xt91228659710.ps tmp/89xt91228659710.png")
> system("convert tmp/9lc0g1228659710.ps tmp/9lc0g1228659710.png")
> system("convert tmp/105lhm1228659710.ps tmp/105lhm1228659710.png")
>
>
> proc.time()
user system elapsed
5.904 3.137 6.313