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,27,0,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,0,51,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 27 0 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 0 51 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
166885 27255 -702 -1348 -11267 -19866
M3 M4 M5 M6 M7 M8
-22365 -26708 -29401 -34672 -37302 -42375
M9 M10 M11
-40194 -6555 1314
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-8108.1 -2398.7 931.7 2110.1 8614.9
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 166884.83 2205.56 75.665 < 2e-16 ***
d 27255.19 5407.21 5.041 7.67e-06 ***
t1 -702.06 35.42 -19.822 < 2e-16 ***
t2 -1348.49 131.30 -10.270 1.73e-13 ***
M1 -11267.27 2565.32 -4.392 6.54e-05 ***
M2 -19865.66 2703.22 -7.349 2.74e-09 ***
M3 -22365.43 2703.02 -8.274 1.17e-10 ***
M4 -26708.44 2711.01 -9.852 6.53e-13 ***
M5 -29401.01 2702.28 -10.880 2.60e-14 ***
M6 -34671.78 2694.69 -12.867 < 2e-16 ***
M7 -37302.35 2688.25 -13.876 < 2e-16 ***
M8 -42374.72 2682.97 -15.794 < 2e-16 ***
M9 -40194.29 2678.86 -15.004 < 2e-16 ***
M10 -6555.46 2675.92 -2.450 0.0182 *
M11 1314.17 2674.15 0.491 0.6255
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4227 on 46 degrees of freedom
Multiple R-squared: 0.9625, Adjusted R-squared: 0.9511
F-statistic: 84.3 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.55192536 0.89614928 0.44807464
[2,] 0.42559277 0.85118555 0.57440723
[3,] 0.33481007 0.66962013 0.66518993
[4,] 0.33749756 0.67499512 0.66250244
[5,] 0.28523131 0.57046262 0.71476869
[6,] 0.22728001 0.45456003 0.77271999
[7,] 0.14938763 0.29877527 0.85061237
[8,] 0.11344011 0.22688022 0.88655989
[9,] 0.07979023 0.15958047 0.92020977
[10,] 0.06530809 0.13061619 0.93469191
[11,] 0.06196378 0.12392757 0.93803622
[12,] 0.05407401 0.10814801 0.94592599
[13,] 0.05290039 0.10580078 0.94709961
[14,] 0.10564176 0.21128352 0.89435824
[15,] 0.07919163 0.15838326 0.92080837
[16,] 0.05119756 0.10239512 0.94880244
[17,] 0.08997203 0.17994405 0.91002797
[18,] 0.14051984 0.28103968 0.85948016
[19,] 0.62591037 0.74817926 0.37408963
[20,] 0.82455724 0.35088552 0.17544276
[21,] 0.96447671 0.07104658 0.03552329
[22,] 0.95255080 0.09489840 0.04744920
[23,] 0.91839788 0.16320425 0.08160212
[24,] 0.84203746 0.31592508 0.15796254
[25,] 0.72469182 0.55061636 0.27530818
[26,] 0.57179071 0.85641858 0.42820929
> postscript(file="/var/www/html/rcomp/tmp/1mnp81229962831.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/2uthw1229962831.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/3sm9q1229962831.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/4pzwm1229962831.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/5y4oj1229962831.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
-7147.5039 -8108.0513 -5494.2213 -1217.1513 -972.5213 -2446.6913 -5021.0613
8 9 10 11 12 13 14
-4735.6313 -4179.0013 -505.7713 1284.6587 2481.8887 2979.2195 2017.6720
15 16 17 18 19 20 21
599.5020 1378.5721 1062.2021 2825.0321 1108.6621 1666.0921 931.7221
22 23 24 25 26 27 28
3844.9521 5178.3821 7185.6120 8614.9428 7758.3954 6547.2254 -4347.9933
29 30 31 32 33 34 35
-2916.9383 -2398.6833 -1315.6283 2123.2268 2647.2818 1828.9368 1688.7918
36 37 38 39 40 41 42
3552.4468 1360.2026 -348.9198 -3569.6648 3073.8303 1820.8853 1777.1403
43 44 45 46 47 48 49
2949.1953 2110.0503 -212.8947 -1425.2397 -1701.3847 -5466.7297 -1825.9739
50 51 52 53 54 55 56
-1319.0963 1917.1587 1112.7422 1006.3722 243.2022 2278.8322 -1163.7378
57 58 59 60 61
812.8922 -3742.8778 -6450.4478 -7753.2178 -3980.8871
> postscript(file="/var/www/html/rcomp/tmp/6ewnc1229962831.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 -7147.5039 NA
1 -8108.0513 -7147.5039
2 -5494.2213 -8108.0513
3 -1217.1513 -5494.2213
4 -972.5213 -1217.1513
5 -2446.6913 -972.5213
6 -5021.0613 -2446.6913
7 -4735.6313 -5021.0613
8 -4179.0013 -4735.6313
9 -505.7713 -4179.0013
10 1284.6587 -505.7713
11 2481.8887 1284.6587
12 2979.2195 2481.8887
13 2017.6720 2979.2195
14 599.5020 2017.6720
15 1378.5721 599.5020
16 1062.2021 1378.5721
17 2825.0321 1062.2021
18 1108.6621 2825.0321
19 1666.0921 1108.6621
20 931.7221 1666.0921
21 3844.9521 931.7221
22 5178.3821 3844.9521
23 7185.6120 5178.3821
24 8614.9428 7185.6120
25 7758.3954 8614.9428
26 6547.2254 7758.3954
27 -4347.9933 6547.2254
28 -2916.9383 -4347.9933
29 -2398.6833 -2916.9383
30 -1315.6283 -2398.6833
31 2123.2268 -1315.6283
32 2647.2818 2123.2268
33 1828.9368 2647.2818
34 1688.7918 1828.9368
35 3552.4468 1688.7918
36 1360.2026 3552.4468
37 -348.9198 1360.2026
38 -3569.6648 -348.9198
39 3073.8303 -3569.6648
40 1820.8853 3073.8303
41 1777.1403 1820.8853
42 2949.1953 1777.1403
43 2110.0503 2949.1953
44 -212.8947 2110.0503
45 -1425.2397 -212.8947
46 -1701.3847 -1425.2397
47 -5466.7297 -1701.3847
48 -1825.9739 -5466.7297
49 -1319.0963 -1825.9739
50 1917.1587 -1319.0963
51 1112.7422 1917.1587
52 1006.3722 1112.7422
53 243.2022 1006.3722
54 2278.8322 243.2022
55 -1163.7378 2278.8322
56 812.8922 -1163.7378
57 -3742.8778 812.8922
58 -6450.4478 -3742.8778
59 -7753.2178 -6450.4478
60 -3980.8871 -7753.2178
61 NA -3980.8871
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -8108.0513 -7147.5039
[2,] -5494.2213 -8108.0513
[3,] -1217.1513 -5494.2213
[4,] -972.5213 -1217.1513
[5,] -2446.6913 -972.5213
[6,] -5021.0613 -2446.6913
[7,] -4735.6313 -5021.0613
[8,] -4179.0013 -4735.6313
[9,] -505.7713 -4179.0013
[10,] 1284.6587 -505.7713
[11,] 2481.8887 1284.6587
[12,] 2979.2195 2481.8887
[13,] 2017.6720 2979.2195
[14,] 599.5020 2017.6720
[15,] 1378.5721 599.5020
[16,] 1062.2021 1378.5721
[17,] 2825.0321 1062.2021
[18,] 1108.6621 2825.0321
[19,] 1666.0921 1108.6621
[20,] 931.7221 1666.0921
[21,] 3844.9521 931.7221
[22,] 5178.3821 3844.9521
[23,] 7185.6120 5178.3821
[24,] 8614.9428 7185.6120
[25,] 7758.3954 8614.9428
[26,] 6547.2254 7758.3954
[27,] -4347.9933 6547.2254
[28,] -2916.9383 -4347.9933
[29,] -2398.6833 -2916.9383
[30,] -1315.6283 -2398.6833
[31,] 2123.2268 -1315.6283
[32,] 2647.2818 2123.2268
[33,] 1828.9368 2647.2818
[34,] 1688.7918 1828.9368
[35,] 3552.4468 1688.7918
[36,] 1360.2026 3552.4468
[37,] -348.9198 1360.2026
[38,] -3569.6648 -348.9198
[39,] 3073.8303 -3569.6648
[40,] 1820.8853 3073.8303
[41,] 1777.1403 1820.8853
[42,] 2949.1953 1777.1403
[43,] 2110.0503 2949.1953
[44,] -212.8947 2110.0503
[45,] -1425.2397 -212.8947
[46,] -1701.3847 -1425.2397
[47,] -5466.7297 -1701.3847
[48,] -1825.9739 -5466.7297
[49,] -1319.0963 -1825.9739
[50,] 1917.1587 -1319.0963
[51,] 1112.7422 1917.1587
[52,] 1006.3722 1112.7422
[53,] 243.2022 1006.3722
[54,] 2278.8322 243.2022
[55,] -1163.7378 2278.8322
[56,] 812.8922 -1163.7378
[57,] -3742.8778 812.8922
[58,] -6450.4478 -3742.8778
[59,] -7753.2178 -6450.4478
[60,] -3980.8871 -7753.2178
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -8108.0513 -7147.5039
2 -5494.2213 -8108.0513
3 -1217.1513 -5494.2213
4 -972.5213 -1217.1513
5 -2446.6913 -972.5213
6 -5021.0613 -2446.6913
7 -4735.6313 -5021.0613
8 -4179.0013 -4735.6313
9 -505.7713 -4179.0013
10 1284.6587 -505.7713
11 2481.8887 1284.6587
12 2979.2195 2481.8887
13 2017.6720 2979.2195
14 599.5020 2017.6720
15 1378.5721 599.5020
16 1062.2021 1378.5721
17 2825.0321 1062.2021
18 1108.6621 2825.0321
19 1666.0921 1108.6621
20 931.7221 1666.0921
21 3844.9521 931.7221
22 5178.3821 3844.9521
23 7185.6120 5178.3821
24 8614.9428 7185.6120
25 7758.3954 8614.9428
26 6547.2254 7758.3954
27 -4347.9933 6547.2254
28 -2916.9383 -4347.9933
29 -2398.6833 -2916.9383
30 -1315.6283 -2398.6833
31 2123.2268 -1315.6283
32 2647.2818 2123.2268
33 1828.9368 2647.2818
34 1688.7918 1828.9368
35 3552.4468 1688.7918
36 1360.2026 3552.4468
37 -348.9198 1360.2026
38 -3569.6648 -348.9198
39 3073.8303 -3569.6648
40 1820.8853 3073.8303
41 1777.1403 1820.8853
42 2949.1953 1777.1403
43 2110.0503 2949.1953
44 -212.8947 2110.0503
45 -1425.2397 -212.8947
46 -1701.3847 -1425.2397
47 -5466.7297 -1701.3847
48 -1825.9739 -5466.7297
49 -1319.0963 -1825.9739
50 1917.1587 -1319.0963
51 1112.7422 1917.1587
52 1006.3722 1112.7422
53 243.2022 1006.3722
54 2278.8322 243.2022
55 -1163.7378 2278.8322
56 812.8922 -1163.7378
57 -3742.8778 812.8922
58 -6450.4478 -3742.8778
59 -7753.2178 -6450.4478
60 -3980.8871 -7753.2178
> 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/7s2qy1229962831.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/8b5jd1229962831.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/9kur61229962831.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/10jmu41229962831.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/11u8d11229962831.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/127ec41229962831.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/13cpfq1229962831.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/14vnem1229962831.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/15yd9b1229962831.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/1605un1229962831.tab")
+ }
> system("convert tmp/1mnp81229962831.ps tmp/1mnp81229962831.png")
> system("convert tmp/2uthw1229962831.ps tmp/2uthw1229962831.png")
> system("convert tmp/3sm9q1229962831.ps tmp/3sm9q1229962831.png")
> system("convert tmp/4pzwm1229962831.ps tmp/4pzwm1229962831.png")
> system("convert tmp/5y4oj1229962831.ps tmp/5y4oj1229962831.png")
> system("convert tmp/6ewnc1229962831.ps tmp/6ewnc1229962831.png")
> system("convert tmp/7s2qy1229962831.ps tmp/7s2qy1229962831.png")
> system("convert tmp/8b5jd1229962831.ps tmp/8b5jd1229962831.png")
> system("convert tmp/9kur61229962831.ps tmp/9kur61229962831.png")
> system("convert tmp/10jmu41229962831.ps tmp/10jmu41229962831.png")
>
>
> proc.time()
user system elapsed
2.492 1.581 3.100