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(313737,312276,309391,302950,300316,304035,333476,337698,335932,323931,313927,314485,313218,309664,302963,298989,298423,301631,329765,335083,327616,309119,295916,291413,291542,284678,276475,272566,264981,263290,296806,303598,286994,276427,266424,267153,268381,262522,255542,253158,243803,250741,280445,285257,270976,261076,255603,260376,263903,264291,263276,262572,256167,264221,293860,300713,287224,275902,271115,277509,279681),dim=c(1,61),dimnames=list(c('HPC'),1:61))
> y <- array(NA,dim=c(1,61),dimnames=list(c('HPC'),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 = '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
HPC M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 313737 1 0 0 0 0 0 0 0 0 0 0 1
2 312276 0 1 0 0 0 0 0 0 0 0 0 2
3 309391 0 0 1 0 0 0 0 0 0 0 0 3
4 302950 0 0 0 1 0 0 0 0 0 0 0 4
5 300316 0 0 0 0 1 0 0 0 0 0 0 5
6 304035 0 0 0 0 0 1 0 0 0 0 0 6
7 333476 0 0 0 0 0 0 1 0 0 0 0 7
8 337698 0 0 0 0 0 0 0 1 0 0 0 8
9 335932 0 0 0 0 0 0 0 0 1 0 0 9
10 323931 0 0 0 0 0 0 0 0 0 1 0 10
11 313927 0 0 0 0 0 0 0 0 0 0 1 11
12 314485 0 0 0 0 0 0 0 0 0 0 0 12
13 313218 1 0 0 0 0 0 0 0 0 0 0 13
14 309664 0 1 0 0 0 0 0 0 0 0 0 14
15 302963 0 0 1 0 0 0 0 0 0 0 0 15
16 298989 0 0 0 1 0 0 0 0 0 0 0 16
17 298423 0 0 0 0 1 0 0 0 0 0 0 17
18 301631 0 0 0 0 0 1 0 0 0 0 0 18
19 329765 0 0 0 0 0 0 1 0 0 0 0 19
20 335083 0 0 0 0 0 0 0 1 0 0 0 20
21 327616 0 0 0 0 0 0 0 0 1 0 0 21
22 309119 0 0 0 0 0 0 0 0 0 1 0 22
23 295916 0 0 0 0 0 0 0 0 0 0 1 23
24 291413 0 0 0 0 0 0 0 0 0 0 0 24
25 291542 1 0 0 0 0 0 0 0 0 0 0 25
26 284678 0 1 0 0 0 0 0 0 0 0 0 26
27 276475 0 0 1 0 0 0 0 0 0 0 0 27
28 272566 0 0 0 1 0 0 0 0 0 0 0 28
29 264981 0 0 0 0 1 0 0 0 0 0 0 29
30 263290 0 0 0 0 0 1 0 0 0 0 0 30
31 296806 0 0 0 0 0 0 1 0 0 0 0 31
32 303598 0 0 0 0 0 0 0 1 0 0 0 32
33 286994 0 0 0 0 0 0 0 0 1 0 0 33
34 276427 0 0 0 0 0 0 0 0 0 1 0 34
35 266424 0 0 0 0 0 0 0 0 0 0 1 35
36 267153 0 0 0 0 0 0 0 0 0 0 0 36
37 268381 1 0 0 0 0 0 0 0 0 0 0 37
38 262522 0 1 0 0 0 0 0 0 0 0 0 38
39 255542 0 0 1 0 0 0 0 0 0 0 0 39
40 253158 0 0 0 1 0 0 0 0 0 0 0 40
41 243803 0 0 0 0 1 0 0 0 0 0 0 41
42 250741 0 0 0 0 0 1 0 0 0 0 0 42
43 280445 0 0 0 0 0 0 1 0 0 0 0 43
44 285257 0 0 0 0 0 0 0 1 0 0 0 44
45 270976 0 0 0 0 0 0 0 0 1 0 0 45
46 261076 0 0 0 0 0 0 0 0 0 1 0 46
47 255603 0 0 0 0 0 0 0 0 0 0 1 47
48 260376 0 0 0 0 0 0 0 0 0 0 0 48
49 263903 1 0 0 0 0 0 0 0 0 0 0 49
50 264291 0 1 0 0 0 0 0 0 0 0 0 50
51 263276 0 0 1 0 0 0 0 0 0 0 0 51
52 262572 0 0 0 1 0 0 0 0 0 0 0 52
53 256167 0 0 0 0 1 0 0 0 0 0 0 53
54 264221 0 0 0 0 0 1 0 0 0 0 0 54
55 293860 0 0 0 0 0 0 1 0 0 0 0 55
56 300713 0 0 0 0 0 0 0 1 0 0 0 56
57 287224 0 0 0 0 0 0 0 0 1 0 0 57
58 275902 0 0 0 0 0 0 0 0 0 1 0 58
59 271115 0 0 0 0 0 0 0 0 0 0 1 59
60 277509 0 0 0 0 0 0 0 0 0 0 0 60
61 279681 1 0 0 0 0 0 0 0 0 0 0 61
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
320655.7 880.3 -6186.7 -10274.9 -12688.7 -16929.2
M6 M7 M8 M9 M10 M11
-11815.0 19340.4 26008.3 16355.5 4966.7 -2658.8
t
-1068.6
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-17950 -11341 960 9075 23328
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 320655.67 6373.96 50.307 < 2e-16 ***
M1 880.29 7433.53 0.118 0.90623
M2 -6186.69 7802.28 -0.793 0.43172
M3 -10274.92 7792.32 -1.319 0.19356
M4 -12688.75 7783.39 -1.630 0.10960
M5 -16929.18 7775.51 -2.177 0.03441 *
M6 -11815.01 7768.67 -1.521 0.13486
M7 19340.36 7762.88 2.491 0.01623 *
M8 26008.33 7758.13 3.352 0.00157 **
M9 16355.49 7754.44 2.109 0.04017 *
M10 4966.66 7751.81 0.641 0.52476
M11 -2658.77 7750.22 -0.343 0.73305
t -1068.57 90.43 -11.816 8.15e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 12250 on 48 degrees of freedom
Multiple R-squared: 0.8044, Adjusted R-squared: 0.7555
F-statistic: 16.45 on 12 and 48 DF, p-value: 4.124e-13
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 2.851324e-03 5.702649e-03 9.971487e-01
[2,] 3.983676e-04 7.967353e-04 9.996016e-01
[3,] 4.971112e-05 9.942224e-05 9.999503e-01
[4,] 6.232425e-06 1.246485e-05 9.999938e-01
[5,] 7.897191e-07 1.579438e-06 9.999992e-01
[6,] 3.360053e-06 6.720106e-06 9.999966e-01
[7,] 1.964309e-04 3.928617e-04 9.998036e-01
[8,] 2.104321e-03 4.208641e-03 9.978957e-01
[9,] 1.535494e-02 3.070988e-02 9.846451e-01
[10,] 2.422183e-02 4.844366e-02 9.757782e-01
[11,] 5.389406e-02 1.077881e-01 9.461059e-01
[12,] 1.068362e-01 2.136724e-01 8.931638e-01
[13,] 1.425096e-01 2.850193e-01 8.574904e-01
[14,] 2.668833e-01 5.337667e-01 7.331167e-01
[15,] 3.914685e-01 7.829370e-01 6.085315e-01
[16,] 4.617088e-01 9.234177e-01 5.382912e-01
[17,] 5.568741e-01 8.862517e-01 4.431259e-01
[18,] 7.646726e-01 4.706548e-01 2.353274e-01
[19,] 9.051121e-01 1.897757e-01 9.488786e-02
[20,] 9.657244e-01 6.855115e-02 3.427558e-02
[21,] 9.871190e-01 2.576191e-02 1.288095e-02
[22,] 9.990632e-01 1.873508e-03 9.367540e-04
[23,] 9.999571e-01 8.571588e-05 4.285794e-05
[24,] 9.999877e-01 2.467108e-05 1.233554e-05
[25,] 9.999980e-01 4.065537e-06 2.032768e-06
[26,] 9.999973e-01 5.421762e-06 2.710881e-06
[27,] 9.999917e-01 1.668011e-05 8.340053e-06
[28,] 9.999924e-01 1.529064e-05 7.645320e-06
[29,] 9.999025e-01 1.949540e-04 9.747698e-05
[30,] 9.988630e-01 2.273989e-03 1.136995e-03
> postscript(file="/var/www/html/rcomp/tmp/11cg71292876106.ps",horizontal=F,onefile=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/21cg71292876106.ps",horizontal=F,onefile=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/34vi51292876107.ps",horizontal=F,onefile=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/44vi51292876107.ps",horizontal=F,onefile=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/54vi51292876107.ps",horizontal=F,onefile=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
-6730.39216 -55.84706 2215.95294 -742.64706 1932.35294 1605.75294
7 8 9 10 11 12
959.95294 -417.44706 8537.95294 8994.35294 7684.35294 6652.15294
13 14 15 16 17 18
5573.43137 10154.97647 8610.77647 8119.17647 12862.17647 12024.57647
19 20 21 22 23 24
10071.77647 9790.37647 13044.77647 7005.17647 2496.17647 -3597.02353
25 26 27 28 29 30
-3279.74510 -2008.20000 -5054.40000 -5481.00000 -7757.00000 -13493.60000
31 32 33 34 35 36
-10064.40000 -8871.80000 -14754.40000 -12864.00000 -14173.00000 -15034.20000
37 38 39 40 41 42
-13617.92157 -11341.37647 -13164.57647 -12066.17647 -16112.17647 -13219.77647
43 44 45 46 47 48
-13602.57647 -14389.97647 -17949.57647 -15392.17647 -12171.17647 -8988.37647
49 50 51 52 53 54
-5273.09804 3250.44706 7392.24706 10170.64706 9074.64706 13083.04706
55 56 57 58 59 60
12635.24706 13888.84706 11121.24706 12256.64706 16163.64706 20967.44706
61
23327.72549
> postscript(file="/var/www/html/rcomp/tmp/6f5h81292876107.ps",horizontal=F,onefile=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 -6730.39216 NA
1 -55.84706 -6730.39216
2 2215.95294 -55.84706
3 -742.64706 2215.95294
4 1932.35294 -742.64706
5 1605.75294 1932.35294
6 959.95294 1605.75294
7 -417.44706 959.95294
8 8537.95294 -417.44706
9 8994.35294 8537.95294
10 7684.35294 8994.35294
11 6652.15294 7684.35294
12 5573.43137 6652.15294
13 10154.97647 5573.43137
14 8610.77647 10154.97647
15 8119.17647 8610.77647
16 12862.17647 8119.17647
17 12024.57647 12862.17647
18 10071.77647 12024.57647
19 9790.37647 10071.77647
20 13044.77647 9790.37647
21 7005.17647 13044.77647
22 2496.17647 7005.17647
23 -3597.02353 2496.17647
24 -3279.74510 -3597.02353
25 -2008.20000 -3279.74510
26 -5054.40000 -2008.20000
27 -5481.00000 -5054.40000
28 -7757.00000 -5481.00000
29 -13493.60000 -7757.00000
30 -10064.40000 -13493.60000
31 -8871.80000 -10064.40000
32 -14754.40000 -8871.80000
33 -12864.00000 -14754.40000
34 -14173.00000 -12864.00000
35 -15034.20000 -14173.00000
36 -13617.92157 -15034.20000
37 -11341.37647 -13617.92157
38 -13164.57647 -11341.37647
39 -12066.17647 -13164.57647
40 -16112.17647 -12066.17647
41 -13219.77647 -16112.17647
42 -13602.57647 -13219.77647
43 -14389.97647 -13602.57647
44 -17949.57647 -14389.97647
45 -15392.17647 -17949.57647
46 -12171.17647 -15392.17647
47 -8988.37647 -12171.17647
48 -5273.09804 -8988.37647
49 3250.44706 -5273.09804
50 7392.24706 3250.44706
51 10170.64706 7392.24706
52 9074.64706 10170.64706
53 13083.04706 9074.64706
54 12635.24706 13083.04706
55 13888.84706 12635.24706
56 11121.24706 13888.84706
57 12256.64706 11121.24706
58 16163.64706 12256.64706
59 20967.44706 16163.64706
60 23327.72549 20967.44706
61 NA 23327.72549
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -55.84706 -6730.39216
[2,] 2215.95294 -55.84706
[3,] -742.64706 2215.95294
[4,] 1932.35294 -742.64706
[5,] 1605.75294 1932.35294
[6,] 959.95294 1605.75294
[7,] -417.44706 959.95294
[8,] 8537.95294 -417.44706
[9,] 8994.35294 8537.95294
[10,] 7684.35294 8994.35294
[11,] 6652.15294 7684.35294
[12,] 5573.43137 6652.15294
[13,] 10154.97647 5573.43137
[14,] 8610.77647 10154.97647
[15,] 8119.17647 8610.77647
[16,] 12862.17647 8119.17647
[17,] 12024.57647 12862.17647
[18,] 10071.77647 12024.57647
[19,] 9790.37647 10071.77647
[20,] 13044.77647 9790.37647
[21,] 7005.17647 13044.77647
[22,] 2496.17647 7005.17647
[23,] -3597.02353 2496.17647
[24,] -3279.74510 -3597.02353
[25,] -2008.20000 -3279.74510
[26,] -5054.40000 -2008.20000
[27,] -5481.00000 -5054.40000
[28,] -7757.00000 -5481.00000
[29,] -13493.60000 -7757.00000
[30,] -10064.40000 -13493.60000
[31,] -8871.80000 -10064.40000
[32,] -14754.40000 -8871.80000
[33,] -12864.00000 -14754.40000
[34,] -14173.00000 -12864.00000
[35,] -15034.20000 -14173.00000
[36,] -13617.92157 -15034.20000
[37,] -11341.37647 -13617.92157
[38,] -13164.57647 -11341.37647
[39,] -12066.17647 -13164.57647
[40,] -16112.17647 -12066.17647
[41,] -13219.77647 -16112.17647
[42,] -13602.57647 -13219.77647
[43,] -14389.97647 -13602.57647
[44,] -17949.57647 -14389.97647
[45,] -15392.17647 -17949.57647
[46,] -12171.17647 -15392.17647
[47,] -8988.37647 -12171.17647
[48,] -5273.09804 -8988.37647
[49,] 3250.44706 -5273.09804
[50,] 7392.24706 3250.44706
[51,] 10170.64706 7392.24706
[52,] 9074.64706 10170.64706
[53,] 13083.04706 9074.64706
[54,] 12635.24706 13083.04706
[55,] 13888.84706 12635.24706
[56,] 11121.24706 13888.84706
[57,] 12256.64706 11121.24706
[58,] 16163.64706 12256.64706
[59,] 20967.44706 16163.64706
[60,] 23327.72549 20967.44706
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -55.84706 -6730.39216
2 2215.95294 -55.84706
3 -742.64706 2215.95294
4 1932.35294 -742.64706
5 1605.75294 1932.35294
6 959.95294 1605.75294
7 -417.44706 959.95294
8 8537.95294 -417.44706
9 8994.35294 8537.95294
10 7684.35294 8994.35294
11 6652.15294 7684.35294
12 5573.43137 6652.15294
13 10154.97647 5573.43137
14 8610.77647 10154.97647
15 8119.17647 8610.77647
16 12862.17647 8119.17647
17 12024.57647 12862.17647
18 10071.77647 12024.57647
19 9790.37647 10071.77647
20 13044.77647 9790.37647
21 7005.17647 13044.77647
22 2496.17647 7005.17647
23 -3597.02353 2496.17647
24 -3279.74510 -3597.02353
25 -2008.20000 -3279.74510
26 -5054.40000 -2008.20000
27 -5481.00000 -5054.40000
28 -7757.00000 -5481.00000
29 -13493.60000 -7757.00000
30 -10064.40000 -13493.60000
31 -8871.80000 -10064.40000
32 -14754.40000 -8871.80000
33 -12864.00000 -14754.40000
34 -14173.00000 -12864.00000
35 -15034.20000 -14173.00000
36 -13617.92157 -15034.20000
37 -11341.37647 -13617.92157
38 -13164.57647 -11341.37647
39 -12066.17647 -13164.57647
40 -16112.17647 -12066.17647
41 -13219.77647 -16112.17647
42 -13602.57647 -13219.77647
43 -14389.97647 -13602.57647
44 -17949.57647 -14389.97647
45 -15392.17647 -17949.57647
46 -12171.17647 -15392.17647
47 -8988.37647 -12171.17647
48 -5273.09804 -8988.37647
49 3250.44706 -5273.09804
50 7392.24706 3250.44706
51 10170.64706 7392.24706
52 9074.64706 10170.64706
53 13083.04706 9074.64706
54 12635.24706 13083.04706
55 13888.84706 12635.24706
56 11121.24706 13888.84706
57 12256.64706 11121.24706
58 16163.64706 12256.64706
59 20967.44706 16163.64706
60 23327.72549 20967.44706
> 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/7peyb1292876107.ps",horizontal=F,onefile=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/8peyb1292876107.ps",horizontal=F,onefile=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/9peyb1292876107.ps",horizontal=F,onefile=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/100nxw1292876107.ps",horizontal=F,onefile=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/11loe21292876107.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/12p6u81292876107.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/13lyaz1292876107.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/14tkyt1292876107.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/15shpa1292876107.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/16or511292876107.tab")
+ }
>
> try(system("convert tmp/11cg71292876106.ps tmp/11cg71292876106.png",intern=TRUE))
character(0)
> try(system("convert tmp/21cg71292876106.ps tmp/21cg71292876106.png",intern=TRUE))
character(0)
> try(system("convert tmp/34vi51292876107.ps tmp/34vi51292876107.png",intern=TRUE))
character(0)
> try(system("convert tmp/44vi51292876107.ps tmp/44vi51292876107.png",intern=TRUE))
character(0)
> try(system("convert tmp/54vi51292876107.ps tmp/54vi51292876107.png",intern=TRUE))
character(0)
> try(system("convert tmp/6f5h81292876107.ps tmp/6f5h81292876107.png",intern=TRUE))
character(0)
> try(system("convert tmp/7peyb1292876107.ps tmp/7peyb1292876107.png",intern=TRUE))
character(0)
> try(system("convert tmp/8peyb1292876107.ps tmp/8peyb1292876107.png",intern=TRUE))
character(0)
> try(system("convert tmp/9peyb1292876107.ps tmp/9peyb1292876107.png",intern=TRUE))
character(0)
> try(system("convert tmp/100nxw1292876107.ps tmp/100nxw1292876107.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.481 1.643 5.949