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(269645,0,267037,0,258113,0,262813,0,267413,0,267366,0,264777,0,258863,0,254844,0,254868,0,277267,0,285351,0,286602,0,283042,0,276687,0,277915,0,277128,0,277103,0,275037,0,270150,0,267140,0,264993,0,287259,0,291186,0,292300,0,288186,0,281477,0,282656,0,280190,0,280408,0,276836,0,275216,0,274352,0,271311,0,289802,0,290726,0,292300,0,278506,0,269826,0,265861,0,269034,0,264176,0,255198,0,253353,0,246057,0,235372,0,258556,0,260993,0,254663,0,250643,0,243422,0,247105,0,248541,0,245039,0,237080,0,237085,0,225554,0,226839,0,247934,0,248333,1,246969,1,245098,1,246263,1,255765,1,264319,1,268347,1,273046,1,273963,1,267430,1,271993,1,292710,1,295881,1),dim=c(2,72),dimnames=list(c('Y','x'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('Y','x'),1:72))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Y x M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 269645 0 1 0 0 0 0 0 0 0 0 0 0
2 267037 0 0 1 0 0 0 0 0 0 0 0 0
3 258113 0 0 0 1 0 0 0 0 0 0 0 0
4 262813 0 0 0 0 1 0 0 0 0 0 0 0
5 267413 0 0 0 0 0 1 0 0 0 0 0 0
6 267366 0 0 0 0 0 0 1 0 0 0 0 0
7 264777 0 0 0 0 0 0 0 1 0 0 0 0
8 258863 0 0 0 0 0 0 0 0 1 0 0 0
9 254844 0 0 0 0 0 0 0 0 0 1 0 0
10 254868 0 0 0 0 0 0 0 0 0 0 1 0
11 277267 0 0 0 0 0 0 0 0 0 0 0 1
12 285351 0 0 0 0 0 0 0 0 0 0 0 0
13 286602 0 1 0 0 0 0 0 0 0 0 0 0
14 283042 0 0 1 0 0 0 0 0 0 0 0 0
15 276687 0 0 0 1 0 0 0 0 0 0 0 0
16 277915 0 0 0 0 1 0 0 0 0 0 0 0
17 277128 0 0 0 0 0 1 0 0 0 0 0 0
18 277103 0 0 0 0 0 0 1 0 0 0 0 0
19 275037 0 0 0 0 0 0 0 1 0 0 0 0
20 270150 0 0 0 0 0 0 0 0 1 0 0 0
21 267140 0 0 0 0 0 0 0 0 0 1 0 0
22 264993 0 0 0 0 0 0 0 0 0 0 1 0
23 287259 0 0 0 0 0 0 0 0 0 0 0 1
24 291186 0 0 0 0 0 0 0 0 0 0 0 0
25 292300 0 1 0 0 0 0 0 0 0 0 0 0
26 288186 0 0 1 0 0 0 0 0 0 0 0 0
27 281477 0 0 0 1 0 0 0 0 0 0 0 0
28 282656 0 0 0 0 1 0 0 0 0 0 0 0
29 280190 0 0 0 0 0 1 0 0 0 0 0 0
30 280408 0 0 0 0 0 0 1 0 0 0 0 0
31 276836 0 0 0 0 0 0 0 1 0 0 0 0
32 275216 0 0 0 0 0 0 0 0 1 0 0 0
33 274352 0 0 0 0 0 0 0 0 0 1 0 0
34 271311 0 0 0 0 0 0 0 0 0 0 1 0
35 289802 0 0 0 0 0 0 0 0 0 0 0 1
36 290726 0 0 0 0 0 0 0 0 0 0 0 0
37 292300 0 1 0 0 0 0 0 0 0 0 0 0
38 278506 0 0 1 0 0 0 0 0 0 0 0 0
39 269826 0 0 0 1 0 0 0 0 0 0 0 0
40 265861 0 0 0 0 1 0 0 0 0 0 0 0
41 269034 0 0 0 0 0 1 0 0 0 0 0 0
42 264176 0 0 0 0 0 0 1 0 0 0 0 0
43 255198 0 0 0 0 0 0 0 1 0 0 0 0
44 253353 0 0 0 0 0 0 0 0 1 0 0 0
45 246057 0 0 0 0 0 0 0 0 0 1 0 0
46 235372 0 0 0 0 0 0 0 0 0 0 1 0
47 258556 0 0 0 0 0 0 0 0 0 0 0 1
48 260993 0 0 0 0 0 0 0 0 0 0 0 0
49 254663 0 1 0 0 0 0 0 0 0 0 0 0
50 250643 0 0 1 0 0 0 0 0 0 0 0 0
51 243422 0 0 0 1 0 0 0 0 0 0 0 0
52 247105 0 0 0 0 1 0 0 0 0 0 0 0
53 248541 0 0 0 0 0 1 0 0 0 0 0 0
54 245039 0 0 0 0 0 0 1 0 0 0 0 0
55 237080 0 0 0 0 0 0 0 1 0 0 0 0
56 237085 0 0 0 0 0 0 0 0 1 0 0 0
57 225554 0 0 0 0 0 0 0 0 0 1 0 0
58 226839 0 0 0 0 0 0 0 0 0 0 1 0
59 247934 0 0 0 0 0 0 0 0 0 0 0 1
60 248333 1 0 0 0 0 0 0 0 0 0 0 0
61 246969 1 1 0 0 0 0 0 0 0 0 0 0
62 245098 1 0 1 0 0 0 0 0 0 0 0 0
63 246263 1 0 0 1 0 0 0 0 0 0 0 0
64 255765 1 0 0 0 1 0 0 0 0 0 0 0
65 264319 1 0 0 0 0 1 0 0 0 0 0 0
66 268347 1 0 0 0 0 0 1 0 0 0 0 0
67 273046 1 0 0 0 0 0 0 1 0 0 0 0
68 273963 1 0 0 0 0 0 0 0 1 0 0 0
69 267430 1 0 0 0 0 0 0 0 0 1 0 0
70 271993 1 0 0 0 0 0 0 0 0 0 1 0
71 292710 1 0 0 0 0 0 0 0 0 0 0 1
72 295881 1 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) x M1 M2 M3 M4
279491 -2239 -5372 -10366 -16487 -13766
M5 M6 M7 M8 M9 M10
-11347 -12045 -15456 -17680 -23222 -24889
M11
-3530
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-30715 -11285 1098 12562 19630
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 279492 6995 39.956 <2e-16 ***
x -2239 5128 -0.437 0.6639
M1 -5372 9631 -0.558 0.5791
M2 -10366 9631 -1.076 0.2861
M3 -16487 9631 -1.712 0.0922 .
M4 -13766 9631 -1.429 0.1582
M5 -11347 9631 -1.178 0.2434
M6 -12045 9631 -1.251 0.2160
M7 -15456 9631 -1.605 0.1139
M8 -17680 9631 -1.836 0.0714 .
M9 -23222 9631 -2.411 0.0190 *
M10 -24889 9631 -2.584 0.0123 *
M11 -3530 9631 -0.367 0.7153
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 16610 on 59 degrees of freedom
Multiple R-squared: 0.1842, Adjusted R-squared: 0.0183
F-statistic: 1.11 on 12 and 59 DF, p-value: 0.3694
> 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.3294000699 0.658800140 0.6705999
[2,] 0.2039696723 0.407939345 0.7960303
[3,] 0.1228367784 0.245673557 0.8771632
[4,] 0.0739785015 0.147957003 0.9260215
[5,] 0.0450370868 0.090074174 0.9549629
[6,] 0.0287509894 0.057501979 0.9712490
[7,] 0.0166961912 0.033392382 0.9833038
[8,] 0.0096289888 0.019257978 0.9903710
[9,] 0.0048723891 0.009744778 0.9951276
[10,] 0.0045882782 0.009176556 0.9954117
[11,] 0.0042541880 0.008508376 0.9957458
[12,] 0.0042247475 0.008449495 0.9957753
[13,] 0.0037164593 0.007432919 0.9962835
[14,] 0.0023941103 0.004788221 0.9976059
[15,] 0.0016218333 0.003243667 0.9983782
[16,] 0.0010752792 0.002150558 0.9989247
[17,] 0.0008478795 0.001695759 0.9991521
[18,] 0.0010076284 0.002015257 0.9989924
[19,] 0.0010691233 0.002138247 0.9989309
[20,] 0.0009010058 0.001802012 0.9990990
[21,] 0.0007861093 0.001572219 0.9992139
[22,] 0.0024828188 0.004965638 0.9975172
[23,] 0.0042266369 0.008453274 0.9957734
[24,] 0.0064523091 0.012904618 0.9935477
[25,] 0.0074277882 0.014855576 0.9925722
[26,] 0.0076541746 0.015308349 0.9923458
[27,] 0.0079409332 0.015881866 0.9920591
[28,] 0.0093492414 0.018698483 0.9906508
[29,] 0.0086379495 0.017275899 0.9913621
[30,] 0.0106567679 0.021313536 0.9893432
[31,] 0.0196128469 0.039225694 0.9803872
[32,] 0.0249271007 0.049854201 0.9750729
[33,] 0.0341501842 0.068300368 0.9658498
[34,] 0.0870241829 0.174048366 0.9129758
[35,] 0.1788006050 0.357601210 0.8211994
[36,] 0.2656481673 0.531296335 0.7343518
[37,] 0.3288405844 0.657681169 0.6711594
[38,] 0.3532082085 0.706416417 0.6467918
[39,] 0.3364532243 0.672906449 0.6635468
[40,] 0.2637327189 0.527465438 0.7362673
[41,] 0.1785784118 0.357156824 0.8214216
> postscript(file="/var/www/html/rcomp/tmp/111m41258740371.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/218r21258740371.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/34rso1258740371.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/45b5u1258740371.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/5p9vl1258740371.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 = 72
Frequency = 1
1 2 3 4 5 6
-4474.73016 -2088.23016 -4891.56349 -2912.73016 -731.06349 -80.39683
7 8 9 10 11 12
741.43651 -2948.56349 -1425.39683 265.43651 1305.76984 5859.53968
13 14 15 16 17 18
12482.26984 13916.76984 13682.43651 12189.26984 8983.93651 9656.60317
19 20 21 22 23 24
11001.43651 8338.43651 10870.60317 10390.43651 11297.76984 11694.53968
25 26 27 28 29 30
18180.26984 19060.76984 18472.43651 16930.26984 12045.93651 12961.60317
31 32 33 34 35 36
12800.43651 13404.43651 18082.60317 16708.43651 13840.76984 11234.53968
37 38 39 40 41 42
18180.26984 9380.76984 6821.43651 135.26984 889.93651 -3270.39683
43 44 45 46 47 48
-8837.56349 -8458.56349 -10212.39683 -19230.56349 -17405.23016 -18498.46032
49 50 51 52 53 54
-19456.73016 -18482.23016 -19582.56349 -18620.73016 -19603.06349 -22407.39683
55 56 57 58 59 60
-26955.56349 -24726.56349 -30715.39683 -27763.56349 -28027.23016 -28919.07937
61 62 63 64 65 66
-24911.34921 -21787.84921 -14502.18254 -7721.34921 -1585.68254 3139.98413
67 68 69 70 71 72
11249.81746 14390.81746 13399.98413 19629.81746 18988.15079 18628.92063
> postscript(file="/var/www/html/rcomp/tmp/6c4v11258740371.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 -4474.73016 NA
1 -2088.23016 -4474.73016
2 -4891.56349 -2088.23016
3 -2912.73016 -4891.56349
4 -731.06349 -2912.73016
5 -80.39683 -731.06349
6 741.43651 -80.39683
7 -2948.56349 741.43651
8 -1425.39683 -2948.56349
9 265.43651 -1425.39683
10 1305.76984 265.43651
11 5859.53968 1305.76984
12 12482.26984 5859.53968
13 13916.76984 12482.26984
14 13682.43651 13916.76984
15 12189.26984 13682.43651
16 8983.93651 12189.26984
17 9656.60317 8983.93651
18 11001.43651 9656.60317
19 8338.43651 11001.43651
20 10870.60317 8338.43651
21 10390.43651 10870.60317
22 11297.76984 10390.43651
23 11694.53968 11297.76984
24 18180.26984 11694.53968
25 19060.76984 18180.26984
26 18472.43651 19060.76984
27 16930.26984 18472.43651
28 12045.93651 16930.26984
29 12961.60317 12045.93651
30 12800.43651 12961.60317
31 13404.43651 12800.43651
32 18082.60317 13404.43651
33 16708.43651 18082.60317
34 13840.76984 16708.43651
35 11234.53968 13840.76984
36 18180.26984 11234.53968
37 9380.76984 18180.26984
38 6821.43651 9380.76984
39 135.26984 6821.43651
40 889.93651 135.26984
41 -3270.39683 889.93651
42 -8837.56349 -3270.39683
43 -8458.56349 -8837.56349
44 -10212.39683 -8458.56349
45 -19230.56349 -10212.39683
46 -17405.23016 -19230.56349
47 -18498.46032 -17405.23016
48 -19456.73016 -18498.46032
49 -18482.23016 -19456.73016
50 -19582.56349 -18482.23016
51 -18620.73016 -19582.56349
52 -19603.06349 -18620.73016
53 -22407.39683 -19603.06349
54 -26955.56349 -22407.39683
55 -24726.56349 -26955.56349
56 -30715.39683 -24726.56349
57 -27763.56349 -30715.39683
58 -28027.23016 -27763.56349
59 -28919.07937 -28027.23016
60 -24911.34921 -28919.07937
61 -21787.84921 -24911.34921
62 -14502.18254 -21787.84921
63 -7721.34921 -14502.18254
64 -1585.68254 -7721.34921
65 3139.98413 -1585.68254
66 11249.81746 3139.98413
67 14390.81746 11249.81746
68 13399.98413 14390.81746
69 19629.81746 13399.98413
70 18988.15079 19629.81746
71 18628.92063 18988.15079
72 NA 18628.92063
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2088.23016 -4474.73016
[2,] -4891.56349 -2088.23016
[3,] -2912.73016 -4891.56349
[4,] -731.06349 -2912.73016
[5,] -80.39683 -731.06349
[6,] 741.43651 -80.39683
[7,] -2948.56349 741.43651
[8,] -1425.39683 -2948.56349
[9,] 265.43651 -1425.39683
[10,] 1305.76984 265.43651
[11,] 5859.53968 1305.76984
[12,] 12482.26984 5859.53968
[13,] 13916.76984 12482.26984
[14,] 13682.43651 13916.76984
[15,] 12189.26984 13682.43651
[16,] 8983.93651 12189.26984
[17,] 9656.60317 8983.93651
[18,] 11001.43651 9656.60317
[19,] 8338.43651 11001.43651
[20,] 10870.60317 8338.43651
[21,] 10390.43651 10870.60317
[22,] 11297.76984 10390.43651
[23,] 11694.53968 11297.76984
[24,] 18180.26984 11694.53968
[25,] 19060.76984 18180.26984
[26,] 18472.43651 19060.76984
[27,] 16930.26984 18472.43651
[28,] 12045.93651 16930.26984
[29,] 12961.60317 12045.93651
[30,] 12800.43651 12961.60317
[31,] 13404.43651 12800.43651
[32,] 18082.60317 13404.43651
[33,] 16708.43651 18082.60317
[34,] 13840.76984 16708.43651
[35,] 11234.53968 13840.76984
[36,] 18180.26984 11234.53968
[37,] 9380.76984 18180.26984
[38,] 6821.43651 9380.76984
[39,] 135.26984 6821.43651
[40,] 889.93651 135.26984
[41,] -3270.39683 889.93651
[42,] -8837.56349 -3270.39683
[43,] -8458.56349 -8837.56349
[44,] -10212.39683 -8458.56349
[45,] -19230.56349 -10212.39683
[46,] -17405.23016 -19230.56349
[47,] -18498.46032 -17405.23016
[48,] -19456.73016 -18498.46032
[49,] -18482.23016 -19456.73016
[50,] -19582.56349 -18482.23016
[51,] -18620.73016 -19582.56349
[52,] -19603.06349 -18620.73016
[53,] -22407.39683 -19603.06349
[54,] -26955.56349 -22407.39683
[55,] -24726.56349 -26955.56349
[56,] -30715.39683 -24726.56349
[57,] -27763.56349 -30715.39683
[58,] -28027.23016 -27763.56349
[59,] -28919.07937 -28027.23016
[60,] -24911.34921 -28919.07937
[61,] -21787.84921 -24911.34921
[62,] -14502.18254 -21787.84921
[63,] -7721.34921 -14502.18254
[64,] -1585.68254 -7721.34921
[65,] 3139.98413 -1585.68254
[66,] 11249.81746 3139.98413
[67,] 14390.81746 11249.81746
[68,] 13399.98413 14390.81746
[69,] 19629.81746 13399.98413
[70,] 18988.15079 19629.81746
[71,] 18628.92063 18988.15079
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2088.23016 -4474.73016
2 -4891.56349 -2088.23016
3 -2912.73016 -4891.56349
4 -731.06349 -2912.73016
5 -80.39683 -731.06349
6 741.43651 -80.39683
7 -2948.56349 741.43651
8 -1425.39683 -2948.56349
9 265.43651 -1425.39683
10 1305.76984 265.43651
11 5859.53968 1305.76984
12 12482.26984 5859.53968
13 13916.76984 12482.26984
14 13682.43651 13916.76984
15 12189.26984 13682.43651
16 8983.93651 12189.26984
17 9656.60317 8983.93651
18 11001.43651 9656.60317
19 8338.43651 11001.43651
20 10870.60317 8338.43651
21 10390.43651 10870.60317
22 11297.76984 10390.43651
23 11694.53968 11297.76984
24 18180.26984 11694.53968
25 19060.76984 18180.26984
26 18472.43651 19060.76984
27 16930.26984 18472.43651
28 12045.93651 16930.26984
29 12961.60317 12045.93651
30 12800.43651 12961.60317
31 13404.43651 12800.43651
32 18082.60317 13404.43651
33 16708.43651 18082.60317
34 13840.76984 16708.43651
35 11234.53968 13840.76984
36 18180.26984 11234.53968
37 9380.76984 18180.26984
38 6821.43651 9380.76984
39 135.26984 6821.43651
40 889.93651 135.26984
41 -3270.39683 889.93651
42 -8837.56349 -3270.39683
43 -8458.56349 -8837.56349
44 -10212.39683 -8458.56349
45 -19230.56349 -10212.39683
46 -17405.23016 -19230.56349
47 -18498.46032 -17405.23016
48 -19456.73016 -18498.46032
49 -18482.23016 -19456.73016
50 -19582.56349 -18482.23016
51 -18620.73016 -19582.56349
52 -19603.06349 -18620.73016
53 -22407.39683 -19603.06349
54 -26955.56349 -22407.39683
55 -24726.56349 -26955.56349
56 -30715.39683 -24726.56349
57 -27763.56349 -30715.39683
58 -28027.23016 -27763.56349
59 -28919.07937 -28027.23016
60 -24911.34921 -28919.07937
61 -21787.84921 -24911.34921
62 -14502.18254 -21787.84921
63 -7721.34921 -14502.18254
64 -1585.68254 -7721.34921
65 3139.98413 -1585.68254
66 11249.81746 3139.98413
67 14390.81746 11249.81746
68 13399.98413 14390.81746
69 19629.81746 13399.98413
70 18988.15079 19629.81746
71 18628.92063 18988.15079
> 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/7bjvc1258740371.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/8dzu81258740371.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/9xm5c1258740371.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/10nwf31258740371.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/116a8q1258740371.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/12h7jr1258740371.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/13kips1258740371.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/14b0iw1258740371.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/15ah801258740372.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/167fjf1258740372.tab")
+ }
>
> system("convert tmp/111m41258740371.ps tmp/111m41258740371.png")
> system("convert tmp/218r21258740371.ps tmp/218r21258740371.png")
> system("convert tmp/34rso1258740371.ps tmp/34rso1258740371.png")
> system("convert tmp/45b5u1258740371.ps tmp/45b5u1258740371.png")
> system("convert tmp/5p9vl1258740371.ps tmp/5p9vl1258740371.png")
> system("convert tmp/6c4v11258740371.ps tmp/6c4v11258740371.png")
> system("convert tmp/7bjvc1258740371.ps tmp/7bjvc1258740371.png")
> system("convert tmp/8dzu81258740371.ps tmp/8dzu81258740371.png")
> system("convert tmp/9xm5c1258740371.ps tmp/9xm5c1258740371.png")
> system("convert tmp/10nwf31258740371.ps tmp/10nwf31258740371.png")
>
>
> proc.time()
user system elapsed
2.552 1.566 2.949