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(115.6,0,120.3,0,121.9,0,121.7,0,118.9,0,113.4,0,114,0,117.5,0,120.9,0,125.1,0,124.7,0,128.2,0,149.7,0,163.6,0,173.9,0,164.5,0,154.2,0,147.9,0,159.3,0,170.3,0,170,0,174.2,0,190.8,0,179.9,0,240.8,0,241.9,0,241.1,0,239.6,0,220.8,0,209.3,0,209.9,0,228.3,0,242.1,0,226.4,0,231.5,0,229.7,0,257.6,0,260,0,264.4,0,268.8,0,271.4,0,273.8,0,277.4,0,268.2,0,264.6,0,266.6,0,266,0,267.4,0,289.8,0,294,0,310.3,0,311.7,0,302.1,0,298.2,0,299.2,0,296.2,0,299,0,300,0,299.4,0,300.2,0,470.2,0,472.1,0,484.8,0,513.4,1,547.2,1,548.1,1,544.7,1,521.1,1,459,1,413.2,1),dim=c(2,70),dimnames=list(c('y','D'),1:70))
> y <- array(NA,dim=c(2,70),dimnames=list(c('y','D'),1:70))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> 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 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 115.6 0 1 0 0 0 0 0 0 0 0 0 0 1
2 120.3 0 0 1 0 0 0 0 0 0 0 0 0 2
3 121.9 0 0 0 1 0 0 0 0 0 0 0 0 3
4 121.7 0 0 0 0 1 0 0 0 0 0 0 0 4
5 118.9 0 0 0 0 0 1 0 0 0 0 0 0 5
6 113.4 0 0 0 0 0 0 1 0 0 0 0 0 6
7 114.0 0 0 0 0 0 0 0 1 0 0 0 0 7
8 117.5 0 0 0 0 0 0 0 0 1 0 0 0 8
9 120.9 0 0 0 0 0 0 0 0 0 1 0 0 9
10 125.1 0 0 0 0 0 0 0 0 0 0 1 0 10
11 124.7 0 0 0 0 0 0 0 0 0 0 0 1 11
12 128.2 0 0 0 0 0 0 0 0 0 0 0 0 12
13 149.7 0 1 0 0 0 0 0 0 0 0 0 0 13
14 163.6 0 0 1 0 0 0 0 0 0 0 0 0 14
15 173.9 0 0 0 1 0 0 0 0 0 0 0 0 15
16 164.5 0 0 0 0 1 0 0 0 0 0 0 0 16
17 154.2 0 0 0 0 0 1 0 0 0 0 0 0 17
18 147.9 0 0 0 0 0 0 1 0 0 0 0 0 18
19 159.3 0 0 0 0 0 0 0 1 0 0 0 0 19
20 170.3 0 0 0 0 0 0 0 0 1 0 0 0 20
21 170.0 0 0 0 0 0 0 0 0 0 1 0 0 21
22 174.2 0 0 0 0 0 0 0 0 0 0 1 0 22
23 190.8 0 0 0 0 0 0 0 0 0 0 0 1 23
24 179.9 0 0 0 0 0 0 0 0 0 0 0 0 24
25 240.8 0 1 0 0 0 0 0 0 0 0 0 0 25
26 241.9 0 0 1 0 0 0 0 0 0 0 0 0 26
27 241.1 0 0 0 1 0 0 0 0 0 0 0 0 27
28 239.6 0 0 0 0 1 0 0 0 0 0 0 0 28
29 220.8 0 0 0 0 0 1 0 0 0 0 0 0 29
30 209.3 0 0 0 0 0 0 1 0 0 0 0 0 30
31 209.9 0 0 0 0 0 0 0 1 0 0 0 0 31
32 228.3 0 0 0 0 0 0 0 0 1 0 0 0 32
33 242.1 0 0 0 0 0 0 0 0 0 1 0 0 33
34 226.4 0 0 0 0 0 0 0 0 0 0 1 0 34
35 231.5 0 0 0 0 0 0 0 0 0 0 0 1 35
36 229.7 0 0 0 0 0 0 0 0 0 0 0 0 36
37 257.6 0 1 0 0 0 0 0 0 0 0 0 0 37
38 260.0 0 0 1 0 0 0 0 0 0 0 0 0 38
39 264.4 0 0 0 1 0 0 0 0 0 0 0 0 39
40 268.8 0 0 0 0 1 0 0 0 0 0 0 0 40
41 271.4 0 0 0 0 0 1 0 0 0 0 0 0 41
42 273.8 0 0 0 0 0 0 1 0 0 0 0 0 42
43 277.4 0 0 0 0 0 0 0 1 0 0 0 0 43
44 268.2 0 0 0 0 0 0 0 0 1 0 0 0 44
45 264.6 0 0 0 0 0 0 0 0 0 1 0 0 45
46 266.6 0 0 0 0 0 0 0 0 0 0 1 0 46
47 266.0 0 0 0 0 0 0 0 0 0 0 0 1 47
48 267.4 0 0 0 0 0 0 0 0 0 0 0 0 48
49 289.8 0 1 0 0 0 0 0 0 0 0 0 0 49
50 294.0 0 0 1 0 0 0 0 0 0 0 0 0 50
51 310.3 0 0 0 1 0 0 0 0 0 0 0 0 51
52 311.7 0 0 0 0 1 0 0 0 0 0 0 0 52
53 302.1 0 0 0 0 0 1 0 0 0 0 0 0 53
54 298.2 0 0 0 0 0 0 1 0 0 0 0 0 54
55 299.2 0 0 0 0 0 0 0 1 0 0 0 0 55
56 296.2 0 0 0 0 0 0 0 0 1 0 0 0 56
57 299.0 0 0 0 0 0 0 0 0 0 1 0 0 57
58 300.0 0 0 0 0 0 0 0 0 0 0 1 0 58
59 299.4 0 0 0 0 0 0 0 0 0 0 0 1 59
60 300.2 0 0 0 0 0 0 0 0 0 0 0 0 60
61 470.2 0 1 0 0 0 0 0 0 0 0 0 0 61
62 472.1 0 0 1 0 0 0 0 0 0 0 0 0 62
63 484.8 0 0 0 1 0 0 0 0 0 0 0 0 63
64 513.4 1 0 0 0 1 0 0 0 0 0 0 0 64
65 547.2 1 0 0 0 0 1 0 0 0 0 0 0 65
66 548.1 1 0 0 0 0 0 1 0 0 0 0 0 66
67 544.7 1 0 0 0 0 0 0 1 0 0 0 0 67
68 521.1 1 0 0 0 0 0 0 0 1 0 0 0 68
69 459.0 1 0 0 0 0 0 0 0 0 1 0 0 69
70 413.2 1 0 0 0 0 0 0 0 0 0 1 0 70
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) D M1 M2 M3 M4
63.838 133.844 54.709 55.041 58.090 35.298
M5 M6 M7 M8 M9 M10
30.081 21.729 19.662 14.810 2.776 -9.942
M11 t
5.768 4.368
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-80.28810 -15.85721 0.06999 11.50494 87.69825
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 63.8379 15.7491 4.053 0.000157 ***
D 133.8436 14.8696 9.001 1.80e-12 ***
M1 54.7092 18.6211 2.938 0.004789 **
M2 55.0413 18.6100 2.958 0.004534 **
M3 58.0902 18.6013 3.123 0.002833 **
M4 35.2984 18.7889 1.879 0.065497 .
M5 30.0806 18.7706 1.603 0.114663
M6 21.7294 18.7547 1.159 0.251531
M7 19.6616 18.7412 1.049 0.298635
M8 14.8104 18.7302 0.791 0.432441
M9 2.7759 18.7216 0.148 0.882661
M10 -9.9419 18.7154 -0.531 0.597370
M11 5.7678 19.4180 0.297 0.767539
t 4.3678 0.2143 20.380 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 30.7 on 56 degrees of freedom
Multiple R-squared: 0.9447, Adjusted R-squared: 0.9319
F-statistic: 73.66 on 13 and 56 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,] 6.310047e-03 1.262009e-02 0.9936900
[2,] 1.208581e-03 2.417163e-03 0.9987914
[3,] 1.932906e-04 3.865812e-04 0.9998067
[4,] 7.075217e-05 1.415043e-04 0.9999292
[5,] 1.342870e-05 2.685739e-05 0.9999866
[6,] 2.524639e-06 5.049277e-06 0.9999975
[7,] 6.880017e-06 1.376003e-05 0.9999931
[8,] 1.431275e-06 2.862550e-06 0.9999986
[9,] 4.042848e-05 8.085697e-05 0.9999596
[10,] 3.140338e-05 6.280677e-05 0.9999686
[11,] 1.137865e-05 2.275731e-05 0.9999886
[12,] 4.852384e-06 9.704769e-06 0.9999951
[13,] 1.238116e-06 2.476232e-06 0.9999988
[14,] 3.295391e-07 6.590783e-07 0.9999997
[15,] 1.040554e-07 2.081109e-07 0.9999999
[16,] 2.523670e-08 5.047340e-08 1.0000000
[17,] 1.617558e-08 3.235115e-08 1.0000000
[18,] 9.027086e-09 1.805417e-08 1.0000000
[19,] 4.065594e-09 8.131187e-09 1.0000000
[20,] 1.794512e-09 3.589024e-09 1.0000000
[21,] 1.561191e-09 3.122383e-09 1.0000000
[22,] 1.563553e-09 3.127106e-09 1.0000000
[23,] 1.168696e-09 2.337393e-09 1.0000000
[24,] 3.676966e-10 7.353932e-10 1.0000000
[25,] 8.568811e-11 1.713762e-10 1.0000000
[26,] 3.209426e-11 6.418852e-11 1.0000000
[27,] 1.126582e-11 2.253164e-11 1.0000000
[28,] 3.552533e-12 7.105066e-12 1.0000000
[29,] 5.901171e-12 1.180234e-11 1.0000000
[30,] 2.041468e-10 4.082937e-10 1.0000000
[31,] 2.966824e-09 5.933649e-09 1.0000000
[32,] 3.089296e-06 6.178592e-06 0.9999969
[33,] 3.559883e-06 7.119767e-06 0.9999964
[34,] 3.192643e-06 6.385287e-06 0.9999968
[35,] 1.092795e-06 2.185590e-06 0.9999989
[36,] 2.965857e-07 5.931714e-07 0.9999997
[37,] 2.667267e-07 5.334534e-07 0.9999997
> postscript(file="/var/www/html/rcomp/tmp/1xsvm1229365078.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/2r9n21229365078.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/3a4zu1229365078.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/4l8ur1229365078.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/5mab81229365078.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 = 70
Frequency = 1
1 2 3 4 5 6
-7.31491228 -7.31491228 -13.13157895 5.09235589 3.14235589 1.62568922
7 8 9 10 11 12
-0.07431078 3.90902256 14.97568922 27.52568922 7.04807018 11.94807018
13 14 15 16 17 18
-25.62894737 -16.42894737 -13.54561404 -4.52167920 -13.97167920 -16.28834586
19 20 21 22 23 24
-7.18834586 4.29498747 11.66165414 24.21165414 20.73403509 11.23403509
25 26 27 28 29 30
13.05701754 9.45701754 1.24035088 18.16428571 0.21428571 -7.30238095
31 32 33 34 35 36
-9.00238095 9.88095238 31.34761905 23.99761905 9.02000000 8.62000000
37 38 39 40 41 42
-22.55701754 -24.85701754 -27.87368421 -5.04974937 -1.59974937 4.78358396
43 44 45 46 47 48
6.08358396 -2.63308271 1.43358396 11.78358396 -8.89403509 -6.09403509
49 50 51 52 53 54
-42.77105263 -43.27105263 -34.38771930 -14.56378446 -23.31378446 -23.23045113
55 56 57 58 59 60
-24.53045113 -27.04711779 -16.58045113 -7.23045113 -27.90807018 -25.70807018
61 62 63 64 65 66
85.21491228 82.41491228 87.69824561 0.87857143 35.52857143 40.41190476
67 68 69 70
34.71190476 11.59523810 -42.83809524 -80.28809524
> postscript(file="/var/www/html/rcomp/tmp/60c0e1229365078.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 = 70
Frequency = 1
lag(myerror, k = 1) myerror
0 -7.31491228 NA
1 -7.31491228 -7.31491228
2 -13.13157895 -7.31491228
3 5.09235589 -13.13157895
4 3.14235589 5.09235589
5 1.62568922 3.14235589
6 -0.07431078 1.62568922
7 3.90902256 -0.07431078
8 14.97568922 3.90902256
9 27.52568922 14.97568922
10 7.04807018 27.52568922
11 11.94807018 7.04807018
12 -25.62894737 11.94807018
13 -16.42894737 -25.62894737
14 -13.54561404 -16.42894737
15 -4.52167920 -13.54561404
16 -13.97167920 -4.52167920
17 -16.28834586 -13.97167920
18 -7.18834586 -16.28834586
19 4.29498747 -7.18834586
20 11.66165414 4.29498747
21 24.21165414 11.66165414
22 20.73403509 24.21165414
23 11.23403509 20.73403509
24 13.05701754 11.23403509
25 9.45701754 13.05701754
26 1.24035088 9.45701754
27 18.16428571 1.24035088
28 0.21428571 18.16428571
29 -7.30238095 0.21428571
30 -9.00238095 -7.30238095
31 9.88095238 -9.00238095
32 31.34761905 9.88095238
33 23.99761905 31.34761905
34 9.02000000 23.99761905
35 8.62000000 9.02000000
36 -22.55701754 8.62000000
37 -24.85701754 -22.55701754
38 -27.87368421 -24.85701754
39 -5.04974937 -27.87368421
40 -1.59974937 -5.04974937
41 4.78358396 -1.59974937
42 6.08358396 4.78358396
43 -2.63308271 6.08358396
44 1.43358396 -2.63308271
45 11.78358396 1.43358396
46 -8.89403509 11.78358396
47 -6.09403509 -8.89403509
48 -42.77105263 -6.09403509
49 -43.27105263 -42.77105263
50 -34.38771930 -43.27105263
51 -14.56378446 -34.38771930
52 -23.31378446 -14.56378446
53 -23.23045113 -23.31378446
54 -24.53045113 -23.23045113
55 -27.04711779 -24.53045113
56 -16.58045113 -27.04711779
57 -7.23045113 -16.58045113
58 -27.90807018 -7.23045113
59 -25.70807018 -27.90807018
60 85.21491228 -25.70807018
61 82.41491228 85.21491228
62 87.69824561 82.41491228
63 0.87857143 87.69824561
64 35.52857143 0.87857143
65 40.41190476 35.52857143
66 34.71190476 40.41190476
67 11.59523810 34.71190476
68 -42.83809524 11.59523810
69 -80.28809524 -42.83809524
70 NA -80.28809524
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -7.31491228 -7.31491228
[2,] -13.13157895 -7.31491228
[3,] 5.09235589 -13.13157895
[4,] 3.14235589 5.09235589
[5,] 1.62568922 3.14235589
[6,] -0.07431078 1.62568922
[7,] 3.90902256 -0.07431078
[8,] 14.97568922 3.90902256
[9,] 27.52568922 14.97568922
[10,] 7.04807018 27.52568922
[11,] 11.94807018 7.04807018
[12,] -25.62894737 11.94807018
[13,] -16.42894737 -25.62894737
[14,] -13.54561404 -16.42894737
[15,] -4.52167920 -13.54561404
[16,] -13.97167920 -4.52167920
[17,] -16.28834586 -13.97167920
[18,] -7.18834586 -16.28834586
[19,] 4.29498747 -7.18834586
[20,] 11.66165414 4.29498747
[21,] 24.21165414 11.66165414
[22,] 20.73403509 24.21165414
[23,] 11.23403509 20.73403509
[24,] 13.05701754 11.23403509
[25,] 9.45701754 13.05701754
[26,] 1.24035088 9.45701754
[27,] 18.16428571 1.24035088
[28,] 0.21428571 18.16428571
[29,] -7.30238095 0.21428571
[30,] -9.00238095 -7.30238095
[31,] 9.88095238 -9.00238095
[32,] 31.34761905 9.88095238
[33,] 23.99761905 31.34761905
[34,] 9.02000000 23.99761905
[35,] 8.62000000 9.02000000
[36,] -22.55701754 8.62000000
[37,] -24.85701754 -22.55701754
[38,] -27.87368421 -24.85701754
[39,] -5.04974937 -27.87368421
[40,] -1.59974937 -5.04974937
[41,] 4.78358396 -1.59974937
[42,] 6.08358396 4.78358396
[43,] -2.63308271 6.08358396
[44,] 1.43358396 -2.63308271
[45,] 11.78358396 1.43358396
[46,] -8.89403509 11.78358396
[47,] -6.09403509 -8.89403509
[48,] -42.77105263 -6.09403509
[49,] -43.27105263 -42.77105263
[50,] -34.38771930 -43.27105263
[51,] -14.56378446 -34.38771930
[52,] -23.31378446 -14.56378446
[53,] -23.23045113 -23.31378446
[54,] -24.53045113 -23.23045113
[55,] -27.04711779 -24.53045113
[56,] -16.58045113 -27.04711779
[57,] -7.23045113 -16.58045113
[58,] -27.90807018 -7.23045113
[59,] -25.70807018 -27.90807018
[60,] 85.21491228 -25.70807018
[61,] 82.41491228 85.21491228
[62,] 87.69824561 82.41491228
[63,] 0.87857143 87.69824561
[64,] 35.52857143 0.87857143
[65,] 40.41190476 35.52857143
[66,] 34.71190476 40.41190476
[67,] 11.59523810 34.71190476
[68,] -42.83809524 11.59523810
[69,] -80.28809524 -42.83809524
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -7.31491228 -7.31491228
2 -13.13157895 -7.31491228
3 5.09235589 -13.13157895
4 3.14235589 5.09235589
5 1.62568922 3.14235589
6 -0.07431078 1.62568922
7 3.90902256 -0.07431078
8 14.97568922 3.90902256
9 27.52568922 14.97568922
10 7.04807018 27.52568922
11 11.94807018 7.04807018
12 -25.62894737 11.94807018
13 -16.42894737 -25.62894737
14 -13.54561404 -16.42894737
15 -4.52167920 -13.54561404
16 -13.97167920 -4.52167920
17 -16.28834586 -13.97167920
18 -7.18834586 -16.28834586
19 4.29498747 -7.18834586
20 11.66165414 4.29498747
21 24.21165414 11.66165414
22 20.73403509 24.21165414
23 11.23403509 20.73403509
24 13.05701754 11.23403509
25 9.45701754 13.05701754
26 1.24035088 9.45701754
27 18.16428571 1.24035088
28 0.21428571 18.16428571
29 -7.30238095 0.21428571
30 -9.00238095 -7.30238095
31 9.88095238 -9.00238095
32 31.34761905 9.88095238
33 23.99761905 31.34761905
34 9.02000000 23.99761905
35 8.62000000 9.02000000
36 -22.55701754 8.62000000
37 -24.85701754 -22.55701754
38 -27.87368421 -24.85701754
39 -5.04974937 -27.87368421
40 -1.59974937 -5.04974937
41 4.78358396 -1.59974937
42 6.08358396 4.78358396
43 -2.63308271 6.08358396
44 1.43358396 -2.63308271
45 11.78358396 1.43358396
46 -8.89403509 11.78358396
47 -6.09403509 -8.89403509
48 -42.77105263 -6.09403509
49 -43.27105263 -42.77105263
50 -34.38771930 -43.27105263
51 -14.56378446 -34.38771930
52 -23.31378446 -14.56378446
53 -23.23045113 -23.31378446
54 -24.53045113 -23.23045113
55 -27.04711779 -24.53045113
56 -16.58045113 -27.04711779
57 -7.23045113 -16.58045113
58 -27.90807018 -7.23045113
59 -25.70807018 -27.90807018
60 85.21491228 -25.70807018
61 82.41491228 85.21491228
62 87.69824561 82.41491228
63 0.87857143 87.69824561
64 35.52857143 0.87857143
65 40.41190476 35.52857143
66 34.71190476 40.41190476
67 11.59523810 34.71190476
68 -42.83809524 11.59523810
69 -80.28809524 -42.83809524
> 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/7sy5f1229365078.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/8ln021229365078.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/9oinx1229365078.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/101ra41229365078.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/11yp2q1229365078.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/1233bq1229365078.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/13rlq31229365078.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/14zxrg1229365078.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/15lz7a1229365078.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/16tvq11229365078.tab")
+ }
>
> system("convert tmp/1xsvm1229365078.ps tmp/1xsvm1229365078.png")
> system("convert tmp/2r9n21229365078.ps tmp/2r9n21229365078.png")
> system("convert tmp/3a4zu1229365078.ps tmp/3a4zu1229365078.png")
> system("convert tmp/4l8ur1229365078.ps tmp/4l8ur1229365078.png")
> system("convert tmp/5mab81229365078.ps tmp/5mab81229365078.png")
> system("convert tmp/60c0e1229365078.ps tmp/60c0e1229365078.png")
> system("convert tmp/7sy5f1229365078.ps tmp/7sy5f1229365078.png")
> system("convert tmp/8ln021229365078.ps tmp/8ln021229365078.png")
> system("convert tmp/9oinx1229365078.ps tmp/9oinx1229365078.png")
> system("convert tmp/101ra41229365078.ps tmp/101ra41229365078.png")
>
>
> proc.time()
user system elapsed
2.577 1.586 3.226