R version 2.9.0 (2009-04-17)
Copyright (C) 2009 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(286602,326011,283042,328282,276687,317480,277915,317539,277128,313737,277103,312276,275037,309391,270150,302950,267140,300316,264993,304035,287259,333476,291186,337698,292300,335932,288186,323931,281477,313927,282656,314485,280190,313218,280408,309664,276836,302963,275216,298989,274352,298423,271311,301631,289802,329765,290726,335083,292300,327616,278506,309119,269826,295916,265861,291413,269034,291542,264176,284678,255198,276475,253353,272566,246057,264981,235372,263290,258556,296806,260993,303598,254663,286994,250643,276427,243422,266424,247105,267153,248541,268381,245039,262522,237080,255542,237085,253158,225554,243803,226839,250741,247934,280445,248333,285257,246969,270976,245098,261076,246263,255603,255765,260376,264319,263903,268347,264291,273046,263276,273963,262572,267430,256167,271993,264221,292710,293860,295881,300713,293299,287224),dim=c(2,61),dimnames=list(c('Y','X'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Y','X'),1:61))
> 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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 286602 326011 1 0 0 0 0 0 0 0 0 0 0
2 283042 328282 0 1 0 0 0 0 0 0 0 0 0
3 276687 317480 0 0 1 0 0 0 0 0 0 0 0
4 277915 317539 0 0 0 1 0 0 0 0 0 0 0
5 277128 313737 0 0 0 0 1 0 0 0 0 0 0
6 277103 312276 0 0 0 0 0 1 0 0 0 0 0
7 275037 309391 0 0 0 0 0 0 1 0 0 0 0
8 270150 302950 0 0 0 0 0 0 0 1 0 0 0
9 267140 300316 0 0 0 0 0 0 0 0 1 0 0
10 264993 304035 0 0 0 0 0 0 0 0 0 1 0
11 287259 333476 0 0 0 0 0 0 0 0 0 0 1
12 291186 337698 0 0 0 0 0 0 0 0 0 0 0
13 292300 335932 1 0 0 0 0 0 0 0 0 0 0
14 288186 323931 0 1 0 0 0 0 0 0 0 0 0
15 281477 313927 0 0 1 0 0 0 0 0 0 0 0
16 282656 314485 0 0 0 1 0 0 0 0 0 0 0
17 280190 313218 0 0 0 0 1 0 0 0 0 0 0
18 280408 309664 0 0 0 0 0 1 0 0 0 0 0
19 276836 302963 0 0 0 0 0 0 1 0 0 0 0
20 275216 298989 0 0 0 0 0 0 0 1 0 0 0
21 274352 298423 0 0 0 0 0 0 0 0 1 0 0
22 271311 301631 0 0 0 0 0 0 0 0 0 1 0
23 289802 329765 0 0 0 0 0 0 0 0 0 0 1
24 290726 335083 0 0 0 0 0 0 0 0 0 0 0
25 292300 327616 1 0 0 0 0 0 0 0 0 0 0
26 278506 309119 0 1 0 0 0 0 0 0 0 0 0
27 269826 295916 0 0 1 0 0 0 0 0 0 0 0
28 265861 291413 0 0 0 1 0 0 0 0 0 0 0
29 269034 291542 0 0 0 0 1 0 0 0 0 0 0
30 264176 284678 0 0 0 0 0 1 0 0 0 0 0
31 255198 276475 0 0 0 0 0 0 1 0 0 0 0
32 253353 272566 0 0 0 0 0 0 0 1 0 0 0
33 246057 264981 0 0 0 0 0 0 0 0 1 0 0
34 235372 263290 0 0 0 0 0 0 0 0 0 1 0
35 258556 296806 0 0 0 0 0 0 0 0 0 0 1
36 260993 303598 0 0 0 0 0 0 0 0 0 0 0
37 254663 286994 1 0 0 0 0 0 0 0 0 0 0
38 250643 276427 0 1 0 0 0 0 0 0 0 0 0
39 243422 266424 0 0 1 0 0 0 0 0 0 0 0
40 247105 267153 0 0 0 1 0 0 0 0 0 0 0
41 248541 268381 0 0 0 0 1 0 0 0 0 0 0
42 245039 262522 0 0 0 0 0 1 0 0 0 0 0
43 237080 255542 0 0 0 0 0 0 1 0 0 0 0
44 237085 253158 0 0 0 0 0 0 0 1 0 0 0
45 225554 243803 0 0 0 0 0 0 0 0 1 0 0
46 226839 250741 0 0 0 0 0 0 0 0 0 1 0
47 247934 280445 0 0 0 0 0 0 0 0 0 0 1
48 248333 285257 0 0 0 0 0 0 0 0 0 0 0
49 246969 270976 1 0 0 0 0 0 0 0 0 0 0
50 245098 261076 0 1 0 0 0 0 0 0 0 0 0
51 246263 255603 0 0 1 0 0 0 0 0 0 0 0
52 255765 260376 0 0 0 1 0 0 0 0 0 0 0
53 264319 263903 0 0 0 0 1 0 0 0 0 0 0
54 268347 264291 0 0 0 0 0 1 0 0 0 0 0
55 273046 263276 0 0 0 0 0 0 1 0 0 0 0
56 273963 262572 0 0 0 0 0 0 0 1 0 0 0
57 267430 256167 0 0 0 0 0 0 0 0 1 0 0
58 271993 264221 0 0 0 0 0 0 0 0 0 1 0
59 292710 293860 0 0 0 0 0 0 0 0 0 0 1
60 295881 300713 0 0 0 0 0 0 0 0 0 0 0
61 293299 287224 1 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
99313.06 0.57 4071.35 -1088.08 -1006.70 1134.47
M5 M6 M7 M8 M9 M10
3137.56 4287.70 3651.92 4150.92 1330.30 -2980.73
M11
1020.11
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-14059 -5628 -2307 2949 26194
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.931e+04 2.110e+04 4.707 2.17e-05 ***
X 5.700e-01 6.547e-02 8.706 1.93e-11 ***
M1 4.071e+03 7.006e+03 0.581 0.564
M2 -1.088e+03 7.351e+03 -0.148 0.883
M3 -1.007e+03 7.452e+03 -0.135 0.893
M4 1.134e+03 7.448e+03 0.152 0.880
M5 3.138e+03 7.448e+03 0.421 0.675
M6 4.288e+03 7.496e+03 0.572 0.570
M7 3.652e+03 7.579e+03 0.482 0.632
M8 4.151e+03 7.643e+03 0.543 0.590
M9 1.330e+03 7.753e+03 0.172 0.864
M10 -2.981e+03 7.668e+03 -0.389 0.699
M11 1.020e+03 7.313e+03 0.139 0.890
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11550 on 48 degrees of freedom
Multiple R-squared: 0.6764, Adjusted R-squared: 0.5956
F-statistic: 8.363 on 12 and 48 DF, p-value: 3.214e-08
> 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,] 3.615741e-02 7.231482e-02 0.9638426
[2,] 1.024312e-02 2.048624e-02 0.9897569
[3,] 2.836274e-03 5.672548e-03 0.9971637
[4,] 6.256564e-04 1.251313e-03 0.9993743
[5,] 2.015364e-04 4.030727e-04 0.9997985
[6,] 1.112066e-04 2.224133e-04 0.9998888
[7,] 4.408882e-05 8.817765e-05 0.9999559
[8,] 9.599975e-06 1.919995e-05 0.9999904
[9,] 1.993935e-06 3.987870e-06 0.9999980
[10,] 4.206264e-07 8.412529e-07 0.9999996
[11,] 5.728640e-07 1.145728e-06 0.9999994
[12,] 2.253710e-07 4.507419e-07 0.9999998
[13,] 1.093000e-07 2.186001e-07 0.9999999
[14,] 2.376976e-08 4.753951e-08 1.0000000
[15,] 6.956741e-09 1.391348e-08 1.0000000
[16,] 5.199464e-09 1.039893e-08 1.0000000
[17,] 2.559506e-09 5.119013e-09 1.0000000
[18,] 1.816368e-09 3.632736e-09 1.0000000
[19,] 3.850476e-09 7.700951e-09 1.0000000
[20,] 8.483970e-09 1.696794e-08 1.0000000
[21,] 4.011538e-08 8.023077e-08 1.0000000
[22,] 2.504886e-07 5.009772e-07 0.9999997
[23,] 1.030718e-06 2.061436e-06 0.9999990
[24,] 1.300703e-05 2.601407e-05 0.9999870
[25,] 2.354227e-04 4.708454e-04 0.9997646
[26,] 4.184470e-02 8.368939e-02 0.9581553
[27,] 3.269579e-01 6.539158e-01 0.6730421
[28,] 6.509829e-01 6.980342e-01 0.3490171
[29,] 7.853592e-01 4.292816e-01 0.2146408
[30,] 7.128487e-01 5.743026e-01 0.2871513
> postscript(file="/var/www/html/rcomp/tmp/1gaxg1258964758.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/216ia1258964758.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/3h60l1258964758.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/40x4o1258964758.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/5v9at1258964758.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 = 61
Frequency = 1
1 2 3 4 5 6
-2611.7593 -2306.8193 -2585.9607 -3532.7642 -4155.6787 -4498.0277
7 8 9 10 11 12
-4283.7751 -5998.3452 -4686.3204 -4642.1553 -3158.6432 -618.1122
13 14 15 16 17 18
-2568.8230 5317.2918 4229.2828 2949.0446 -797.8438 295.8370
19 20 21 22 23 24
1179.2456 1325.4622 3604.7074 3046.1474 1499.6618 412.4625
25 26 27 28 29 30
2171.3756 4080.2717 2844.7229 -694.6975 401.6809 -1693.9070
31 32 33 34 35 36
-5360.3443 -5476.1782 -5628.0367 -11038.1206 -10959.3970 -11373.7902
37 38 39 40 41 42
-12310.7008 -5147.9796 -6748.5586 -5622.2684 -6889.3304 -8201.7778
43 44 45 46 47 48
-11546.3366 -10681.4350 -14059.3768 -12418.0721 -12255.4725 -13579.2470
49 50 51 52 53 54
-10874.2895 -1942.7646 2260.5136 6900.6856 11441.1719 14097.8755
55 56 57 58 59 60
20011.2104 20830.4961 20769.0265 25052.2006 24873.8509 25158.6870
61
26194.1970
> postscript(file="/var/www/html/rcomp/tmp/6tzeh1258964758.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 -2611.7593 NA
1 -2306.8193 -2611.7593
2 -2585.9607 -2306.8193
3 -3532.7642 -2585.9607
4 -4155.6787 -3532.7642
5 -4498.0277 -4155.6787
6 -4283.7751 -4498.0277
7 -5998.3452 -4283.7751
8 -4686.3204 -5998.3452
9 -4642.1553 -4686.3204
10 -3158.6432 -4642.1553
11 -618.1122 -3158.6432
12 -2568.8230 -618.1122
13 5317.2918 -2568.8230
14 4229.2828 5317.2918
15 2949.0446 4229.2828
16 -797.8438 2949.0446
17 295.8370 -797.8438
18 1179.2456 295.8370
19 1325.4622 1179.2456
20 3604.7074 1325.4622
21 3046.1474 3604.7074
22 1499.6618 3046.1474
23 412.4625 1499.6618
24 2171.3756 412.4625
25 4080.2717 2171.3756
26 2844.7229 4080.2717
27 -694.6975 2844.7229
28 401.6809 -694.6975
29 -1693.9070 401.6809
30 -5360.3443 -1693.9070
31 -5476.1782 -5360.3443
32 -5628.0367 -5476.1782
33 -11038.1206 -5628.0367
34 -10959.3970 -11038.1206
35 -11373.7902 -10959.3970
36 -12310.7008 -11373.7902
37 -5147.9796 -12310.7008
38 -6748.5586 -5147.9796
39 -5622.2684 -6748.5586
40 -6889.3304 -5622.2684
41 -8201.7778 -6889.3304
42 -11546.3366 -8201.7778
43 -10681.4350 -11546.3366
44 -14059.3768 -10681.4350
45 -12418.0721 -14059.3768
46 -12255.4725 -12418.0721
47 -13579.2470 -12255.4725
48 -10874.2895 -13579.2470
49 -1942.7646 -10874.2895
50 2260.5136 -1942.7646
51 6900.6856 2260.5136
52 11441.1719 6900.6856
53 14097.8755 11441.1719
54 20011.2104 14097.8755
55 20830.4961 20011.2104
56 20769.0265 20830.4961
57 25052.2006 20769.0265
58 24873.8509 25052.2006
59 25158.6870 24873.8509
60 26194.1970 25158.6870
61 NA 26194.1970
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2306.8193 -2611.7593
[2,] -2585.9607 -2306.8193
[3,] -3532.7642 -2585.9607
[4,] -4155.6787 -3532.7642
[5,] -4498.0277 -4155.6787
[6,] -4283.7751 -4498.0277
[7,] -5998.3452 -4283.7751
[8,] -4686.3204 -5998.3452
[9,] -4642.1553 -4686.3204
[10,] -3158.6432 -4642.1553
[11,] -618.1122 -3158.6432
[12,] -2568.8230 -618.1122
[13,] 5317.2918 -2568.8230
[14,] 4229.2828 5317.2918
[15,] 2949.0446 4229.2828
[16,] -797.8438 2949.0446
[17,] 295.8370 -797.8438
[18,] 1179.2456 295.8370
[19,] 1325.4622 1179.2456
[20,] 3604.7074 1325.4622
[21,] 3046.1474 3604.7074
[22,] 1499.6618 3046.1474
[23,] 412.4625 1499.6618
[24,] 2171.3756 412.4625
[25,] 4080.2717 2171.3756
[26,] 2844.7229 4080.2717
[27,] -694.6975 2844.7229
[28,] 401.6809 -694.6975
[29,] -1693.9070 401.6809
[30,] -5360.3443 -1693.9070
[31,] -5476.1782 -5360.3443
[32,] -5628.0367 -5476.1782
[33,] -11038.1206 -5628.0367
[34,] -10959.3970 -11038.1206
[35,] -11373.7902 -10959.3970
[36,] -12310.7008 -11373.7902
[37,] -5147.9796 -12310.7008
[38,] -6748.5586 -5147.9796
[39,] -5622.2684 -6748.5586
[40,] -6889.3304 -5622.2684
[41,] -8201.7778 -6889.3304
[42,] -11546.3366 -8201.7778
[43,] -10681.4350 -11546.3366
[44,] -14059.3768 -10681.4350
[45,] -12418.0721 -14059.3768
[46,] -12255.4725 -12418.0721
[47,] -13579.2470 -12255.4725
[48,] -10874.2895 -13579.2470
[49,] -1942.7646 -10874.2895
[50,] 2260.5136 -1942.7646
[51,] 6900.6856 2260.5136
[52,] 11441.1719 6900.6856
[53,] 14097.8755 11441.1719
[54,] 20011.2104 14097.8755
[55,] 20830.4961 20011.2104
[56,] 20769.0265 20830.4961
[57,] 25052.2006 20769.0265
[58,] 24873.8509 25052.2006
[59,] 25158.6870 24873.8509
[60,] 26194.1970 25158.6870
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2306.8193 -2611.7593
2 -2585.9607 -2306.8193
3 -3532.7642 -2585.9607
4 -4155.6787 -3532.7642
5 -4498.0277 -4155.6787
6 -4283.7751 -4498.0277
7 -5998.3452 -4283.7751
8 -4686.3204 -5998.3452
9 -4642.1553 -4686.3204
10 -3158.6432 -4642.1553
11 -618.1122 -3158.6432
12 -2568.8230 -618.1122
13 5317.2918 -2568.8230
14 4229.2828 5317.2918
15 2949.0446 4229.2828
16 -797.8438 2949.0446
17 295.8370 -797.8438
18 1179.2456 295.8370
19 1325.4622 1179.2456
20 3604.7074 1325.4622
21 3046.1474 3604.7074
22 1499.6618 3046.1474
23 412.4625 1499.6618
24 2171.3756 412.4625
25 4080.2717 2171.3756
26 2844.7229 4080.2717
27 -694.6975 2844.7229
28 401.6809 -694.6975
29 -1693.9070 401.6809
30 -5360.3443 -1693.9070
31 -5476.1782 -5360.3443
32 -5628.0367 -5476.1782
33 -11038.1206 -5628.0367
34 -10959.3970 -11038.1206
35 -11373.7902 -10959.3970
36 -12310.7008 -11373.7902
37 -5147.9796 -12310.7008
38 -6748.5586 -5147.9796
39 -5622.2684 -6748.5586
40 -6889.3304 -5622.2684
41 -8201.7778 -6889.3304
42 -11546.3366 -8201.7778
43 -10681.4350 -11546.3366
44 -14059.3768 -10681.4350
45 -12418.0721 -14059.3768
46 -12255.4725 -12418.0721
47 -13579.2470 -12255.4725
48 -10874.2895 -13579.2470
49 -1942.7646 -10874.2895
50 2260.5136 -1942.7646
51 6900.6856 2260.5136
52 11441.1719 6900.6856
53 14097.8755 11441.1719
54 20011.2104 14097.8755
55 20830.4961 20011.2104
56 20769.0265 20830.4961
57 25052.2006 20769.0265
58 24873.8509 25052.2006
59 25158.6870 24873.8509
60 26194.1970 25158.6870
> 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/7p4zb1258964758.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/8iuif1258964758.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/9bloo1258964758.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/106pe21258964758.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/11hlgq1258964758.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/12wpww1258964758.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/1325i21258964758.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/14iauc1258964758.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/15hmcw1258964758.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/16lwmt1258964758.tab")
+ }
>
> system("convert tmp/1gaxg1258964758.ps tmp/1gaxg1258964758.png")
> system("convert tmp/216ia1258964758.ps tmp/216ia1258964758.png")
> system("convert tmp/3h60l1258964758.ps tmp/3h60l1258964758.png")
> system("convert tmp/40x4o1258964758.ps tmp/40x4o1258964758.png")
> system("convert tmp/5v9at1258964758.ps tmp/5v9at1258964758.png")
> system("convert tmp/6tzeh1258964758.ps tmp/6tzeh1258964758.png")
> system("convert tmp/7p4zb1258964758.ps tmp/7p4zb1258964758.png")
> system("convert tmp/8iuif1258964758.ps tmp/8iuif1258964758.png")
> system("convert tmp/9bloo1258964758.ps tmp/9bloo1258964758.png")
> system("convert tmp/106pe21258964758.ps tmp/106pe21258964758.png")
>
>
> proc.time()
user system elapsed
2.365 1.525 12.487