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,1,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,0,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 1 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 0 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
167137.8 27869.9 -710.6 -1368.4 -11297.1 -19942.2
M3 M4 M5 M6 M7 M8
-25586.2 -26813.2 -29492.7 -34750.4 -37367.9 -42427.1
M9 M10 M11
-40233.6 -6581.7 1301.1
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-8267.4 -1688.1 897.8 1961.0 8605.8
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 167137.82 2103.68 79.450 < 2e-16 ***
d 27869.92 5022.08 5.549 1.36e-06 ***
t1 -710.62 32.89 -21.604 < 2e-16 ***
t2 -1368.40 125.37 -10.914 2.34e-14 ***
M1 -11297.09 2451.34 -4.609 3.23e-05 ***
M2 -19942.17 2581.57 -7.725 7.55e-10 ***
M3 -25586.18 2600.62 -9.838 6.81e-13 ***
M4 -26813.25 2591.12 -10.348 1.36e-13 ***
M5 -29492.72 2582.70 -11.419 5.07e-15 ***
M6 -34750.39 2575.39 -13.493 < 2e-16 ***
M7 -37367.86 2569.18 -14.545 < 2e-16 ***
M8 -42427.12 2564.09 -16.547 < 2e-16 ***
M9 -40233.59 2560.13 -15.715 < 2e-16 ***
M10 -6581.66 2557.29 -2.574 0.0133 *
M11 1301.07 2555.59 0.509 0.6131
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4040 on 46 degrees of freedom
Multiple R-squared: 0.9657, Adjusted R-squared: 0.9553
F-statistic: 92.62 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.58546339 0.82907322 0.41453661
[2,] 0.45981601 0.91963202 0.54018399
[3,] 0.36718267 0.73436535 0.63281733
[4,] 0.37477982 0.74955964 0.62522018
[5,] 0.32348069 0.64696138 0.67651931
[6,] 0.26575872 0.53151744 0.73424128
[7,] 0.18066091 0.36132183 0.81933909
[8,] 0.14210743 0.28421485 0.85789257
[9,] 0.12972512 0.25945025 0.87027488
[10,] 0.08186398 0.16372796 0.91813602
[11,] 0.09445097 0.18890193 0.90554903
[12,] 0.12897997 0.25795995 0.87102003
[13,] 0.18194085 0.36388169 0.81805915
[14,] 0.59421320 0.81157360 0.40578680
[15,] 0.66645575 0.66708851 0.33354425
[16,] 0.67432633 0.65134735 0.32567367
[17,] 0.68380899 0.63238202 0.31619101
[18,] 0.69193277 0.61613446 0.30806723
[19,] 0.97793243 0.04413514 0.02206757
[20,] 0.97747019 0.04505962 0.02252981
[21,] 0.96269786 0.07460427 0.03730214
[22,] 0.95051049 0.09897901 0.04948951
[23,] 0.91576575 0.16846850 0.08423425
[24,] 0.83790031 0.32419939 0.16209969
[25,] 0.71915282 0.56169436 0.28084718
[26,] 0.56647746 0.86704508 0.43352254
> postscript(file="/var/www/html/rcomp/tmp/1eidf1229510939.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/2pdpt1229510939.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/30fjw1229510939.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/4znxx1229510939.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/5le111229510939.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
-7362.1084 -8267.4043 -2500.7793 -1331.0907 -1091.0020 -2569.7134 -5148.6248
8 9 10 11 12 13 14
-4867.7361 -4315.6475 -646.9589 1138.9298 2331.6184 2867.3287 1961.0327
15 16 17 18 19 20 21
3695.6577 1367.3464 1046.4350 2804.7236 1083.8123 1636.7009 897.7895
22 23 24 25 26 27 28
3806.4782 5135.3668 7138.0554 8605.7657 7804.4698 -363.8127 -4553.3457
29 30 31 32 33 34 35
-3115.4786 -2590.4116 -1500.5445 1945.1225 2475.9896 1664.4566 1531.1237
36 37 38 39 40 41 42
3401.5908 1259.0795 -383.4381 -440.0346 3107.4324 1861.2995 1824.3665
43 44 45 46 47 48 49
3003.2336 2170.9006 -145.2323 -1350.7653 -1620.0982 -5378.6312 -1688.1425
50 51 52 53 54 55 56
-1114.6600 -391.0311 1409.6575 1298.7462 531.0348 2562.1234 -884.9879
57 58 59 60 61
1087.1007 -3473.2107 -6185.3220 -7492.6334 -3681.9231
> postscript(file="/var/www/html/rcomp/tmp/6oa7w1229510939.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 -7362.1084 NA
1 -8267.4043 -7362.1084
2 -2500.7793 -8267.4043
3 -1331.0907 -2500.7793
4 -1091.0020 -1331.0907
5 -2569.7134 -1091.0020
6 -5148.6248 -2569.7134
7 -4867.7361 -5148.6248
8 -4315.6475 -4867.7361
9 -646.9589 -4315.6475
10 1138.9298 -646.9589
11 2331.6184 1138.9298
12 2867.3287 2331.6184
13 1961.0327 2867.3287
14 3695.6577 1961.0327
15 1367.3464 3695.6577
16 1046.4350 1367.3464
17 2804.7236 1046.4350
18 1083.8123 2804.7236
19 1636.7009 1083.8123
20 897.7895 1636.7009
21 3806.4782 897.7895
22 5135.3668 3806.4782
23 7138.0554 5135.3668
24 8605.7657 7138.0554
25 7804.4698 8605.7657
26 -363.8127 7804.4698
27 -4553.3457 -363.8127
28 -3115.4786 -4553.3457
29 -2590.4116 -3115.4786
30 -1500.5445 -2590.4116
31 1945.1225 -1500.5445
32 2475.9896 1945.1225
33 1664.4566 2475.9896
34 1531.1237 1664.4566
35 3401.5908 1531.1237
36 1259.0795 3401.5908
37 -383.4381 1259.0795
38 -440.0346 -383.4381
39 3107.4324 -440.0346
40 1861.2995 3107.4324
41 1824.3665 1861.2995
42 3003.2336 1824.3665
43 2170.9006 3003.2336
44 -145.2323 2170.9006
45 -1350.7653 -145.2323
46 -1620.0982 -1350.7653
47 -5378.6312 -1620.0982
48 -1688.1425 -5378.6312
49 -1114.6600 -1688.1425
50 -391.0311 -1114.6600
51 1409.6575 -391.0311
52 1298.7462 1409.6575
53 531.0348 1298.7462
54 2562.1234 531.0348
55 -884.9879 2562.1234
56 1087.1007 -884.9879
57 -3473.2107 1087.1007
58 -6185.3220 -3473.2107
59 -7492.6334 -6185.3220
60 -3681.9231 -7492.6334
61 NA -3681.9231
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -8267.4043 -7362.1084
[2,] -2500.7793 -8267.4043
[3,] -1331.0907 -2500.7793
[4,] -1091.0020 -1331.0907
[5,] -2569.7134 -1091.0020
[6,] -5148.6248 -2569.7134
[7,] -4867.7361 -5148.6248
[8,] -4315.6475 -4867.7361
[9,] -646.9589 -4315.6475
[10,] 1138.9298 -646.9589
[11,] 2331.6184 1138.9298
[12,] 2867.3287 2331.6184
[13,] 1961.0327 2867.3287
[14,] 3695.6577 1961.0327
[15,] 1367.3464 3695.6577
[16,] 1046.4350 1367.3464
[17,] 2804.7236 1046.4350
[18,] 1083.8123 2804.7236
[19,] 1636.7009 1083.8123
[20,] 897.7895 1636.7009
[21,] 3806.4782 897.7895
[22,] 5135.3668 3806.4782
[23,] 7138.0554 5135.3668
[24,] 8605.7657 7138.0554
[25,] 7804.4698 8605.7657
[26,] -363.8127 7804.4698
[27,] -4553.3457 -363.8127
[28,] -3115.4786 -4553.3457
[29,] -2590.4116 -3115.4786
[30,] -1500.5445 -2590.4116
[31,] 1945.1225 -1500.5445
[32,] 2475.9896 1945.1225
[33,] 1664.4566 2475.9896
[34,] 1531.1237 1664.4566
[35,] 3401.5908 1531.1237
[36,] 1259.0795 3401.5908
[37,] -383.4381 1259.0795
[38,] -440.0346 -383.4381
[39,] 3107.4324 -440.0346
[40,] 1861.2995 3107.4324
[41,] 1824.3665 1861.2995
[42,] 3003.2336 1824.3665
[43,] 2170.9006 3003.2336
[44,] -145.2323 2170.9006
[45,] -1350.7653 -145.2323
[46,] -1620.0982 -1350.7653
[47,] -5378.6312 -1620.0982
[48,] -1688.1425 -5378.6312
[49,] -1114.6600 -1688.1425
[50,] -391.0311 -1114.6600
[51,] 1409.6575 -391.0311
[52,] 1298.7462 1409.6575
[53,] 531.0348 1298.7462
[54,] 2562.1234 531.0348
[55,] -884.9879 2562.1234
[56,] 1087.1007 -884.9879
[57,] -3473.2107 1087.1007
[58,] -6185.3220 -3473.2107
[59,] -7492.6334 -6185.3220
[60,] -3681.9231 -7492.6334
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -8267.4043 -7362.1084
2 -2500.7793 -8267.4043
3 -1331.0907 -2500.7793
4 -1091.0020 -1331.0907
5 -2569.7134 -1091.0020
6 -5148.6248 -2569.7134
7 -4867.7361 -5148.6248
8 -4315.6475 -4867.7361
9 -646.9589 -4315.6475
10 1138.9298 -646.9589
11 2331.6184 1138.9298
12 2867.3287 2331.6184
13 1961.0327 2867.3287
14 3695.6577 1961.0327
15 1367.3464 3695.6577
16 1046.4350 1367.3464
17 2804.7236 1046.4350
18 1083.8123 2804.7236
19 1636.7009 1083.8123
20 897.7895 1636.7009
21 3806.4782 897.7895
22 5135.3668 3806.4782
23 7138.0554 5135.3668
24 8605.7657 7138.0554
25 7804.4698 8605.7657
26 -363.8127 7804.4698
27 -4553.3457 -363.8127
28 -3115.4786 -4553.3457
29 -2590.4116 -3115.4786
30 -1500.5445 -2590.4116
31 1945.1225 -1500.5445
32 2475.9896 1945.1225
33 1664.4566 2475.9896
34 1531.1237 1664.4566
35 3401.5908 1531.1237
36 1259.0795 3401.5908
37 -383.4381 1259.0795
38 -440.0346 -383.4381
39 3107.4324 -440.0346
40 1861.2995 3107.4324
41 1824.3665 1861.2995
42 3003.2336 1824.3665
43 2170.9006 3003.2336
44 -145.2323 2170.9006
45 -1350.7653 -145.2323
46 -1620.0982 -1350.7653
47 -5378.6312 -1620.0982
48 -1688.1425 -5378.6312
49 -1114.6600 -1688.1425
50 -391.0311 -1114.6600
51 1409.6575 -391.0311
52 1298.7462 1409.6575
53 531.0348 1298.7462
54 2562.1234 531.0348
55 -884.9879 2562.1234
56 1087.1007 -884.9879
57 -3473.2107 1087.1007
58 -6185.3220 -3473.2107
59 -7492.6334 -6185.3220
60 -3681.9231 -7492.6334
> 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/7prpw1229510939.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/85ycs1229510939.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/9f3ce1229510939.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/10r0kv1229510939.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/111tz51229510939.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/12ae6a1229510939.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/130x4a1229510939.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/14xfbn1229510939.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/15rdry1229510939.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/16t0ms1229510940.tab")
+ }
> system("convert tmp/1eidf1229510939.ps tmp/1eidf1229510939.png")
> system("convert tmp/2pdpt1229510939.ps tmp/2pdpt1229510939.png")
> system("convert tmp/30fjw1229510939.ps tmp/30fjw1229510939.png")
> system("convert tmp/4znxx1229510939.ps tmp/4znxx1229510939.png")
> system("convert tmp/5le111229510939.ps tmp/5le111229510939.png")
> system("convert tmp/6oa7w1229510939.ps tmp/6oa7w1229510939.png")
> system("convert tmp/7prpw1229510939.ps tmp/7prpw1229510939.png")
> system("convert tmp/85ycs1229510939.ps tmp/85ycs1229510939.png")
> system("convert tmp/9f3ce1229510939.ps tmp/9f3ce1229510939.png")
> system("convert tmp/10r0kv1229510939.ps tmp/10r0kv1229510939.png")
>
>
> proc.time()
user system elapsed
2.712 1.760 5.450