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(147768,0,1,0,137507,0,2,0,136919,0,3,0,136151,0,4,0,133001,0,5,0,125554,0,6,0,119647,0,7,0,114158,0,8,0,116193,0,9,0,152803,0,10,0,161761,0,11,0,160942,0,12,0,149470,0,13,0,139208,0,14,0,134588,0,15,0,130322,0,16,0,126611,0,17,0,122401,0,18,0,117352,0,19,0,112135,0,20,0,112879,0,21,0,148729,0,22,0,157230,0,23,0,157221,0,24,0,146681,0,25,0,136524,0,26,0,132111,0,0,27,125326,1,0,28,122716,1,0,29,116615,1,0,30,113719,1,0,31,110737,1,0,32,112093,1,0,33,143565,1,0,34,149946,1,0,35,149147,1,0,36,134339,1,0,37,122683,1,0,38,115614,1,0,39,116566,1,0,40,111272,1,0,41,104609,1,0,42,101802,1,0,43,94542,1,0,44,93051,1,0,45,124129,1,0,46,130374,1,0,47,123946,1,0,48,114971,1,0,49,105531,1,0,50,104919,1,51,0,104782,0,52,0,101281,0,53,0,94545,0,54,0,93248,0,55,0,84031,0,56,0,87486,0,57,0,115867,0,58,0,120327,0,59,0,117008,0,60,0,108811,0,61,0),dim=c(4,61),dimnames=list(c('y','d','t1','t2'),1:61))
> y <- array(NA,dim=c(4,61),dimnames=list(c('y','d','t1','t2'),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 d t1 t2 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 147768 0 1 0 1 0 0 0 0 0 0 0 0 0 0
2 137507 0 2 0 0 1 0 0 0 0 0 0 0 0 0
3 136919 0 3 0 0 0 1 0 0 0 0 0 0 0 0
4 136151 0 4 0 0 0 0 1 0 0 0 0 0 0 0
5 133001 0 5 0 0 0 0 0 1 0 0 0 0 0 0
6 125554 0 6 0 0 0 0 0 0 1 0 0 0 0 0
7 119647 0 7 0 0 0 0 0 0 0 1 0 0 0 0
8 114158 0 8 0 0 0 0 0 0 0 0 1 0 0 0
9 116193 0 9 0 0 0 0 0 0 0 0 0 1 0 0
10 152803 0 10 0 0 0 0 0 0 0 0 0 0 1 0
11 161761 0 11 0 0 0 0 0 0 0 0 0 0 0 1
12 160942 0 12 0 0 0 0 0 0 0 0 0 0 0 0
13 149470 0 13 0 1 0 0 0 0 0 0 0 0 0 0
14 139208 0 14 0 0 1 0 0 0 0 0 0 0 0 0
15 134588 0 15 0 0 0 1 0 0 0 0 0 0 0 0
16 130322 0 16 0 0 0 0 1 0 0 0 0 0 0 0
17 126611 0 17 0 0 0 0 0 1 0 0 0 0 0 0
18 122401 0 18 0 0 0 0 0 0 1 0 0 0 0 0
19 117352 0 19 0 0 0 0 0 0 0 1 0 0 0 0
20 112135 0 20 0 0 0 0 0 0 0 0 1 0 0 0
21 112879 0 21 0 0 0 0 0 0 0 0 0 1 0 0
22 148729 0 22 0 0 0 0 0 0 0 0 0 0 1 0
23 157230 0 23 0 0 0 0 0 0 0 0 0 0 0 1
24 157221 0 24 0 0 0 0 0 0 0 0 0 0 0 0
25 146681 0 25 0 1 0 0 0 0 0 0 0 0 0 0
26 136524 0 26 0 0 1 0 0 0 0 0 0 0 0 0
27 132111 0 0 27 0 0 1 0 0 0 0 0 0 0 0
28 125326 1 0 28 0 0 0 1 0 0 0 0 0 0 0
29 122716 1 0 29 0 0 0 0 1 0 0 0 0 0 0
30 116615 1 0 30 0 0 0 0 0 1 0 0 0 0 0
31 113719 1 0 31 0 0 0 0 0 0 1 0 0 0 0
32 110737 1 0 32 0 0 0 0 0 0 0 1 0 0 0
33 112093 1 0 33 0 0 0 0 0 0 0 0 1 0 0
34 143565 1 0 34 0 0 0 0 0 0 0 0 0 1 0
35 149946 1 0 35 0 0 0 0 0 0 0 0 0 0 1
36 149147 1 0 36 0 0 0 0 0 0 0 0 0 0 0
37 134339 1 0 37 1 0 0 0 0 0 0 0 0 0 0
38 122683 1 0 38 0 1 0 0 0 0 0 0 0 0 0
39 115614 1 0 39 0 0 1 0 0 0 0 0 0 0 0
40 116566 1 0 40 0 0 0 1 0 0 0 0 0 0 0
41 111272 1 0 41 0 0 0 0 1 0 0 0 0 0 0
42 104609 1 0 42 0 0 0 0 0 1 0 0 0 0 0
43 101802 1 0 43 0 0 0 0 0 0 1 0 0 0 0
44 94542 1 0 44 0 0 0 0 0 0 0 1 0 0 0
45 93051 1 0 45 0 0 0 0 0 0 0 0 1 0 0
46 124129 1 0 46 0 0 0 0 0 0 0 0 0 1 0
47 130374 1 0 47 0 0 0 0 0 0 0 0 0 0 1
48 123946 1 0 48 0 0 0 0 0 0 0 0 0 0 0
49 114971 1 0 49 1 0 0 0 0 0 0 0 0 0 0
50 105531 1 0 50 0 1 0 0 0 0 0 0 0 0 0
51 104919 1 51 0 0 0 1 0 0 0 0 0 0 0 0
52 104782 0 52 0 0 0 0 1 0 0 0 0 0 0 0
53 101281 0 53 0 0 0 0 0 1 0 0 0 0 0 0
54 94545 0 54 0 0 0 0 0 0 1 0 0 0 0 0
55 93248 0 55 0 0 0 0 0 0 0 1 0 0 0 0
56 84031 0 56 0 0 0 0 0 0 0 0 1 0 0 0
57 87486 0 57 0 0 0 0 0 0 0 0 0 1 0 0
58 115867 0 58 0 0 0 0 0 0 0 0 0 0 1 0
59 120327 0 59 0 0 0 0 0 0 0 0 0 0 0 1
60 117008 0 60 0 0 0 0 0 0 0 0 0 0 0 0
61 108811 0 61 0 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) d t1 t2 M1 M2
167775.4 5022.9 -747.2 -820.5 -11561.6 -20776.0
M3 M4 M5 M6 M7 M8
-23811.6 -25235.8 -28112.5 -33567.3 -36382.0 -41638.4
M9 M10 M11
-39642.1 -6187.3 1498.2
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-9467 -3400 668 2670 10301
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 167775.45 2663.44 62.992 < 2e-16 ***
d 5022.89 3658.96 1.373 0.176481
t1 -747.24 41.21 -18.135 < 2e-16 ***
t2 -820.52 98.59 -8.322 9.95e-11 ***
M1 -11561.57 3103.76 -3.725 0.000532 ***
M2 -20776.01 3263.35 -6.366 8.16e-08 ***
M3 -23811.59 3271.18 -7.279 3.48e-09 ***
M4 -25235.83 3263.76 -7.732 7.36e-10 ***
M5 -28112.48 3257.21 -8.631 3.54e-11 ***
M6 -33567.32 3251.52 -10.324 1.46e-13 ***
M7 -36381.97 3246.69 -11.206 9.64e-15 ***
M8 -41638.42 3242.74 -12.841 < 2e-16 ***
M9 -39642.06 3239.66 -12.236 4.56e-16 ***
M10 -6187.31 3237.46 -1.911 0.062227 .
M11 1498.25 3236.14 0.463 0.645567
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5116 on 46 degrees of freedom
Multiple R-squared: 0.9451, Adjusted R-squared: 0.9283
F-statistic: 56.51 on 14 and 46 DF, p-value: < 2.2e-16
> 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,] 0.40115691 0.80231382 0.59884309
[2,] 0.29148878 0.58297756 0.70851122
[3,] 0.22666635 0.45333271 0.77333365
[4,] 0.24906709 0.49813418 0.75093291
[5,] 0.23263429 0.46526857 0.76736571
[6,] 0.20930323 0.41860647 0.79069677
[7,] 0.14659019 0.29318038 0.85340981
[8,] 0.12967823 0.25935645 0.87032177
[9,] 0.12064932 0.24129864 0.87935068
[10,] 0.07439174 0.14878349 0.92560826
[11,] 0.09282143 0.18564287 0.90717857
[12,] 0.10751155 0.21502310 0.89248845
[13,] 0.14315585 0.28631171 0.85684415
[14,] 0.49265362 0.98530725 0.50734638
[15,] 0.62768022 0.74463955 0.37231978
[16,] 0.69059580 0.61880839 0.30940420
[17,] 0.63102252 0.73795496 0.36897748
[18,] 0.62365136 0.75269729 0.37634864
[19,] 0.94278071 0.11443858 0.05721929
[20,] 0.95074588 0.09850824 0.04925412
[21,] 0.98482514 0.03034971 0.01517486
[22,] 0.98938554 0.02122892 0.01061446
[23,] 0.98067706 0.03864589 0.01932294
[24,] 0.96149487 0.07701026 0.03850513
[25,] 0.92160869 0.15678262 0.07839131
[26,] 0.82910093 0.34179813 0.17089907
> postscript(file="/var/www/html/rcomp/tmp/1iksz1229962329.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/2z0c91229962329.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/3i46c1229962329.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/4qjc41229962329.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/5pa9x1229962329.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 7
-7698.6318 -7997.9469 -4803.1320 -3399.6424 -2925.7529 -4170.6633 -6515.7738
8 9 10 11 12 13 14
-6001.0842 -5215.1947 -1312.7051 706.9844 2133.4740 2970.2899 2669.9748
15 16 17 18 19 20 21
1832.7898 -261.7207 -348.8311 1643.2584 156.1480 942.8375 437.7271
22 23 24 25 26 27 28
3580.2166 5142.9062 7379.3957 9148.2116 8952.8966 10301.1668 738.0417
29 30 31 32 33 34 35
1825.2074 1999.5730 2738.7387 5833.7044 6013.8700 4851.6357 4367.6014
36 37 38 39 40 41 42
5887.3671 3461.4591 1840.4202 -1372.4888 1824.2769 227.4426 -160.1918
43 44 45 46 47 48 49
667.9739 -515.0604 -3181.8947 -4738.1291 -5358.1634 -9467.3977 -6060.3057
50 51 52 53 54 55 56
-5465.3446 -5958.3358 1099.0445 1221.9341 688.0236 2952.9132 -260.3973
57 58 59 60 61
1945.4923 -2381.0181 -4859.3286 -5932.8390 -1821.0231
> postscript(file="/var/www/html/rcomp/tmp/65zry1229962329.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 -7698.6318 NA
1 -7997.9469 -7698.6318
2 -4803.1320 -7997.9469
3 -3399.6424 -4803.1320
4 -2925.7529 -3399.6424
5 -4170.6633 -2925.7529
6 -6515.7738 -4170.6633
7 -6001.0842 -6515.7738
8 -5215.1947 -6001.0842
9 -1312.7051 -5215.1947
10 706.9844 -1312.7051
11 2133.4740 706.9844
12 2970.2899 2133.4740
13 2669.9748 2970.2899
14 1832.7898 2669.9748
15 -261.7207 1832.7898
16 -348.8311 -261.7207
17 1643.2584 -348.8311
18 156.1480 1643.2584
19 942.8375 156.1480
20 437.7271 942.8375
21 3580.2166 437.7271
22 5142.9062 3580.2166
23 7379.3957 5142.9062
24 9148.2116 7379.3957
25 8952.8966 9148.2116
26 10301.1668 8952.8966
27 738.0417 10301.1668
28 1825.2074 738.0417
29 1999.5730 1825.2074
30 2738.7387 1999.5730
31 5833.7044 2738.7387
32 6013.8700 5833.7044
33 4851.6357 6013.8700
34 4367.6014 4851.6357
35 5887.3671 4367.6014
36 3461.4591 5887.3671
37 1840.4202 3461.4591
38 -1372.4888 1840.4202
39 1824.2769 -1372.4888
40 227.4426 1824.2769
41 -160.1918 227.4426
42 667.9739 -160.1918
43 -515.0604 667.9739
44 -3181.8947 -515.0604
45 -4738.1291 -3181.8947
46 -5358.1634 -4738.1291
47 -9467.3977 -5358.1634
48 -6060.3057 -9467.3977
49 -5465.3446 -6060.3057
50 -5958.3358 -5465.3446
51 1099.0445 -5958.3358
52 1221.9341 1099.0445
53 688.0236 1221.9341
54 2952.9132 688.0236
55 -260.3973 2952.9132
56 1945.4923 -260.3973
57 -2381.0181 1945.4923
58 -4859.3286 -2381.0181
59 -5932.8390 -4859.3286
60 -1821.0231 -5932.8390
61 NA -1821.0231
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -7997.9469 -7698.6318
[2,] -4803.1320 -7997.9469
[3,] -3399.6424 -4803.1320
[4,] -2925.7529 -3399.6424
[5,] -4170.6633 -2925.7529
[6,] -6515.7738 -4170.6633
[7,] -6001.0842 -6515.7738
[8,] -5215.1947 -6001.0842
[9,] -1312.7051 -5215.1947
[10,] 706.9844 -1312.7051
[11,] 2133.4740 706.9844
[12,] 2970.2899 2133.4740
[13,] 2669.9748 2970.2899
[14,] 1832.7898 2669.9748
[15,] -261.7207 1832.7898
[16,] -348.8311 -261.7207
[17,] 1643.2584 -348.8311
[18,] 156.1480 1643.2584
[19,] 942.8375 156.1480
[20,] 437.7271 942.8375
[21,] 3580.2166 437.7271
[22,] 5142.9062 3580.2166
[23,] 7379.3957 5142.9062
[24,] 9148.2116 7379.3957
[25,] 8952.8966 9148.2116
[26,] 10301.1668 8952.8966
[27,] 738.0417 10301.1668
[28,] 1825.2074 738.0417
[29,] 1999.5730 1825.2074
[30,] 2738.7387 1999.5730
[31,] 5833.7044 2738.7387
[32,] 6013.8700 5833.7044
[33,] 4851.6357 6013.8700
[34,] 4367.6014 4851.6357
[35,] 5887.3671 4367.6014
[36,] 3461.4591 5887.3671
[37,] 1840.4202 3461.4591
[38,] -1372.4888 1840.4202
[39,] 1824.2769 -1372.4888
[40,] 227.4426 1824.2769
[41,] -160.1918 227.4426
[42,] 667.9739 -160.1918
[43,] -515.0604 667.9739
[44,] -3181.8947 -515.0604
[45,] -4738.1291 -3181.8947
[46,] -5358.1634 -4738.1291
[47,] -9467.3977 -5358.1634
[48,] -6060.3057 -9467.3977
[49,] -5465.3446 -6060.3057
[50,] -5958.3358 -5465.3446
[51,] 1099.0445 -5958.3358
[52,] 1221.9341 1099.0445
[53,] 688.0236 1221.9341
[54,] 2952.9132 688.0236
[55,] -260.3973 2952.9132
[56,] 1945.4923 -260.3973
[57,] -2381.0181 1945.4923
[58,] -4859.3286 -2381.0181
[59,] -5932.8390 -4859.3286
[60,] -1821.0231 -5932.8390
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -7997.9469 -7698.6318
2 -4803.1320 -7997.9469
3 -3399.6424 -4803.1320
4 -2925.7529 -3399.6424
5 -4170.6633 -2925.7529
6 -6515.7738 -4170.6633
7 -6001.0842 -6515.7738
8 -5215.1947 -6001.0842
9 -1312.7051 -5215.1947
10 706.9844 -1312.7051
11 2133.4740 706.9844
12 2970.2899 2133.4740
13 2669.9748 2970.2899
14 1832.7898 2669.9748
15 -261.7207 1832.7898
16 -348.8311 -261.7207
17 1643.2584 -348.8311
18 156.1480 1643.2584
19 942.8375 156.1480
20 437.7271 942.8375
21 3580.2166 437.7271
22 5142.9062 3580.2166
23 7379.3957 5142.9062
24 9148.2116 7379.3957
25 8952.8966 9148.2116
26 10301.1668 8952.8966
27 738.0417 10301.1668
28 1825.2074 738.0417
29 1999.5730 1825.2074
30 2738.7387 1999.5730
31 5833.7044 2738.7387
32 6013.8700 5833.7044
33 4851.6357 6013.8700
34 4367.6014 4851.6357
35 5887.3671 4367.6014
36 3461.4591 5887.3671
37 1840.4202 3461.4591
38 -1372.4888 1840.4202
39 1824.2769 -1372.4888
40 227.4426 1824.2769
41 -160.1918 227.4426
42 667.9739 -160.1918
43 -515.0604 667.9739
44 -3181.8947 -515.0604
45 -4738.1291 -3181.8947
46 -5358.1634 -4738.1291
47 -9467.3977 -5358.1634
48 -6060.3057 -9467.3977
49 -5465.3446 -6060.3057
50 -5958.3358 -5465.3446
51 1099.0445 -5958.3358
52 1221.9341 1099.0445
53 688.0236 1221.9341
54 2952.9132 688.0236
55 -260.3973 2952.9132
56 1945.4923 -260.3973
57 -2381.0181 1945.4923
58 -4859.3286 -2381.0181
59 -5932.8390 -4859.3286
60 -1821.0231 -5932.8390
> 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/7cof21229962329.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/8gaow1229962329.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/9f1uw1229962329.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/10v8ub1229962329.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/11dwsq1229962329.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/12hqnz1229962329.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/13rjrq1229962329.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/14jie01229962329.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/151bnv1229962329.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/166ni11229962329.tab")
+ }
> system("convert tmp/1iksz1229962329.ps tmp/1iksz1229962329.png")
> system("convert tmp/2z0c91229962329.ps tmp/2z0c91229962329.png")
> system("convert tmp/3i46c1229962329.ps tmp/3i46c1229962329.png")
> system("convert tmp/4qjc41229962329.ps tmp/4qjc41229962329.png")
> system("convert tmp/5pa9x1229962329.ps tmp/5pa9x1229962329.png")
> system("convert tmp/65zry1229962329.ps tmp/65zry1229962329.png")
> system("convert tmp/7cof21229962329.ps tmp/7cof21229962329.png")
> system("convert tmp/8gaow1229962329.ps tmp/8gaow1229962329.png")
> system("convert tmp/9f1uw1229962329.ps tmp/9f1uw1229962329.png")
> system("convert tmp/10v8ub1229962329.ps tmp/10v8ub1229962329.png")
>
>
> proc.time()
user system elapsed
2.518 1.624 3.208