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(31.75,0,27.85,0,27.33,0,29.11,0,28.17,0,28.93,0,32.33,0,34.74,0,33.70,0,34.35,0,35.35,0,34.44,0,33.70,0,32.39,0,28.30,0,29.11,0,28.67,0,28.18,0,29.28,0,29.73,0,26.26,0,26.82,0,27.72,0,27.10,0,27.03,0,25.98,0,25.72,0,25.93,0,24.94,0,21.70,0,17.90,0,17.06,0,16.41,0,16.68,0,18.24,0,16.41,0,15.71,0,13.95,0,12.22,0,14.91,0,14.61,0,15.01,0,15.57,0,16.07,0,15.39,0,15.16,0,15.44,0,15.70,0,17.57,0,18.42,0,17.93,0,18.42,0,17.61,0,17.98,0,17.78,0,17.74,0,19.04,0,19.85,0,20.23,0,20.23,0,21.07,0,21.28,0,21.83,0,21.83,0,22.22,0,22.68,0,23.58,0,23.73,0,23.68,0,23.92,0,24.85,0,26.28,0,27.75,0,29.59,0,29.26,0,29.25,0,28.68,0,26.05,0,27.11,0,29.53,0,31.01,0,32.95,0,32.09,0,31.74,0,32.50,0,33.60,0,32.47,0,34.38,0,32.31,0,30.71,0,30.26,0,27.20,0,24.85,0,22.27,1,18.11,1,18.30,1),dim=c(2,96),dimnames=list(c('x','y'),1:96))
> y <- array(NA,dim=c(2,96),dimnames=list(c('x','y'),1:96))
> 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'
> #'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
x y M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 31.75 0 1 0 0 0 0 0 0 0 0 0 0 1
2 27.85 0 0 1 0 0 0 0 0 0 0 0 0 2
3 27.33 0 0 0 1 0 0 0 0 0 0 0 0 3
4 29.11 0 0 0 0 1 0 0 0 0 0 0 0 4
5 28.17 0 0 0 0 0 1 0 0 0 0 0 0 5
6 28.93 0 0 0 0 0 0 1 0 0 0 0 0 6
7 32.33 0 0 0 0 0 0 0 1 0 0 0 0 7
8 34.74 0 0 0 0 0 0 0 0 1 0 0 0 8
9 33.70 0 0 0 0 0 0 0 0 0 1 0 0 9
10 34.35 0 0 0 0 0 0 0 0 0 0 1 0 10
11 35.35 0 0 0 0 0 0 0 0 0 0 0 1 11
12 34.44 0 0 0 0 0 0 0 0 0 0 0 0 12
13 33.70 0 1 0 0 0 0 0 0 0 0 0 0 13
14 32.39 0 0 1 0 0 0 0 0 0 0 0 0 14
15 28.30 0 0 0 1 0 0 0 0 0 0 0 0 15
16 29.11 0 0 0 0 1 0 0 0 0 0 0 0 16
17 28.67 0 0 0 0 0 1 0 0 0 0 0 0 17
18 28.18 0 0 0 0 0 0 1 0 0 0 0 0 18
19 29.28 0 0 0 0 0 0 0 1 0 0 0 0 19
20 29.73 0 0 0 0 0 0 0 0 1 0 0 0 20
21 26.26 0 0 0 0 0 0 0 0 0 1 0 0 21
22 26.82 0 0 0 0 0 0 0 0 0 0 1 0 22
23 27.72 0 0 0 0 0 0 0 0 0 0 0 1 23
24 27.10 0 0 0 0 0 0 0 0 0 0 0 0 24
25 27.03 0 1 0 0 0 0 0 0 0 0 0 0 25
26 25.98 0 0 1 0 0 0 0 0 0 0 0 0 26
27 25.72 0 0 0 1 0 0 0 0 0 0 0 0 27
28 25.93 0 0 0 0 1 0 0 0 0 0 0 0 28
29 24.94 0 0 0 0 0 1 0 0 0 0 0 0 29
30 21.70 0 0 0 0 0 0 1 0 0 0 0 0 30
31 17.90 0 0 0 0 0 0 0 1 0 0 0 0 31
32 17.06 0 0 0 0 0 0 0 0 1 0 0 0 32
33 16.41 0 0 0 0 0 0 0 0 0 1 0 0 33
34 16.68 0 0 0 0 0 0 0 0 0 0 1 0 34
35 18.24 0 0 0 0 0 0 0 0 0 0 0 1 35
36 16.41 0 0 0 0 0 0 0 0 0 0 0 0 36
37 15.71 0 1 0 0 0 0 0 0 0 0 0 0 37
38 13.95 0 0 1 0 0 0 0 0 0 0 0 0 38
39 12.22 0 0 0 1 0 0 0 0 0 0 0 0 39
40 14.91 0 0 0 0 1 0 0 0 0 0 0 0 40
41 14.61 0 0 0 0 0 1 0 0 0 0 0 0 41
42 15.01 0 0 0 0 0 0 1 0 0 0 0 0 42
43 15.57 0 0 0 0 0 0 0 1 0 0 0 0 43
44 16.07 0 0 0 0 0 0 0 0 1 0 0 0 44
45 15.39 0 0 0 0 0 0 0 0 0 1 0 0 45
46 15.16 0 0 0 0 0 0 0 0 0 0 1 0 46
47 15.44 0 0 0 0 0 0 0 0 0 0 0 1 47
48 15.70 0 0 0 0 0 0 0 0 0 0 0 0 48
49 17.57 0 1 0 0 0 0 0 0 0 0 0 0 49
50 18.42 0 0 1 0 0 0 0 0 0 0 0 0 50
51 17.93 0 0 0 1 0 0 0 0 0 0 0 0 51
52 18.42 0 0 0 0 1 0 0 0 0 0 0 0 52
53 17.61 0 0 0 0 0 1 0 0 0 0 0 0 53
54 17.98 0 0 0 0 0 0 1 0 0 0 0 0 54
55 17.78 0 0 0 0 0 0 0 1 0 0 0 0 55
56 17.74 0 0 0 0 0 0 0 0 1 0 0 0 56
57 19.04 0 0 0 0 0 0 0 0 0 1 0 0 57
58 19.85 0 0 0 0 0 0 0 0 0 0 1 0 58
59 20.23 0 0 0 0 0 0 0 0 0 0 0 1 59
60 20.23 0 0 0 0 0 0 0 0 0 0 0 0 60
61 21.07 0 1 0 0 0 0 0 0 0 0 0 0 61
62 21.28 0 0 1 0 0 0 0 0 0 0 0 0 62
63 21.83 0 0 0 1 0 0 0 0 0 0 0 0 63
64 21.83 0 0 0 0 1 0 0 0 0 0 0 0 64
65 22.22 0 0 0 0 0 1 0 0 0 0 0 0 65
66 22.68 0 0 0 0 0 0 1 0 0 0 0 0 66
67 23.58 0 0 0 0 0 0 0 1 0 0 0 0 67
68 23.73 0 0 0 0 0 0 0 0 1 0 0 0 68
69 23.68 0 0 0 0 0 0 0 0 0 1 0 0 69
70 23.92 0 0 0 0 0 0 0 0 0 0 1 0 70
71 24.85 0 0 0 0 0 0 0 0 0 0 0 1 71
72 26.28 0 0 0 0 0 0 0 0 0 0 0 0 72
73 27.75 0 1 0 0 0 0 0 0 0 0 0 0 73
74 29.59 0 0 1 0 0 0 0 0 0 0 0 0 74
75 29.26 0 0 0 1 0 0 0 0 0 0 0 0 75
76 29.25 0 0 0 0 1 0 0 0 0 0 0 0 76
77 28.68 0 0 0 0 0 1 0 0 0 0 0 0 77
78 26.05 0 0 0 0 0 0 1 0 0 0 0 0 78
79 27.11 0 0 0 0 0 0 0 1 0 0 0 0 79
80 29.53 0 0 0 0 0 0 0 0 1 0 0 0 80
81 31.01 0 0 0 0 0 0 0 0 0 1 0 0 81
82 32.95 0 0 0 0 0 0 0 0 0 0 1 0 82
83 32.09 0 0 0 0 0 0 0 0 0 0 0 1 83
84 31.74 0 0 0 0 0 0 0 0 0 0 0 0 84
85 32.50 0 1 0 0 0 0 0 0 0 0 0 0 85
86 33.60 0 0 1 0 0 0 0 0 0 0 0 0 86
87 32.47 0 0 0 1 0 0 0 0 0 0 0 0 87
88 34.38 0 0 0 0 1 0 0 0 0 0 0 0 88
89 32.31 0 0 0 0 0 1 0 0 0 0 0 0 89
90 30.71 0 0 0 0 0 0 1 0 0 0 0 0 90
91 30.26 0 0 0 0 0 0 0 1 0 0 0 0 91
92 27.20 0 0 0 0 0 0 0 0 1 0 0 0 92
93 24.85 0 0 0 0 0 0 0 0 0 1 0 0 93
94 22.27 1 0 0 0 0 0 0 0 0 0 1 0 94
95 18.11 1 0 0 0 0 0 0 0 0 0 0 1 95
96 18.30 1 0 0 0 0 0 0 0 0 0 0 0 96
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) y M1 M2 M3 M4
24.538057 -4.850693 1.471739 0.972141 -0.024957 0.962946
M5 M6 M7 M8 M9 M10
0.249598 -0.493750 -0.169598 0.082054 -0.597543 0.219196
M11 t
0.225848 -0.002902
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-12.180 -6.389 1.136 4.918 10.618
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 24.538057 2.731479 8.983 7.76e-14 ***
y -4.850693 4.327308 -1.121 0.266
M1 1.471739 3.398493 0.433 0.666
M2 0.972141 3.397596 0.286 0.776
M3 -0.024957 3.396898 -0.007 0.994
M4 0.962946 3.396400 0.284 0.777
M5 0.249598 3.396101 0.073 0.942
M6 -0.493750 3.396001 -0.145 0.885
M7 -0.169598 3.396101 -0.050 0.960
M8 0.082054 3.396400 0.024 0.981
M9 -0.597543 3.396898 -0.176 0.861
M10 0.219196 3.356684 0.065 0.948
M11 0.225848 3.356381 0.067 0.947
t -0.002902 0.026025 -0.112 0.911
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.713 on 82 degrees of freedom
Multiple R-squared: 0.02876, Adjusted R-squared: -0.1252
F-statistic: 0.1868 on 13 and 82 DF, p-value: 0.9991
> 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.02330400 4.660800e-02 9.766960e-01
[2,] 0.00969153 1.938306e-02 9.903085e-01
[3,] 0.01171646 2.343292e-02 9.882835e-01
[4,] 0.02320546 4.641092e-02 9.767945e-01
[5,] 0.05853834 1.170767e-01 9.414617e-01
[6,] 0.08816298 1.763260e-01 9.118370e-01
[7,] 0.12429764 2.485953e-01 8.757024e-01
[8,] 0.16214968 3.242994e-01 8.378503e-01
[9,] 0.16596049 3.319210e-01 8.340395e-01
[10,] 0.16567732 3.313546e-01 8.343227e-01
[11,] 0.20911466 4.182293e-01 7.908853e-01
[12,] 0.27435748 5.487150e-01 7.256425e-01
[13,] 0.41607779 8.321556e-01 5.839222e-01
[14,] 0.61113773 7.777245e-01 3.888623e-01
[15,] 0.90201105 1.959779e-01 9.798895e-02
[16,] 0.99071679 1.856642e-02 9.283210e-03
[17,] 0.99842152 3.156962e-03 1.578481e-03
[18,] 0.99944681 1.106387e-03 5.531936e-04
[19,] 0.99990705 1.858970e-04 9.294849e-05
[20,] 0.99997431 5.137884e-05 2.568942e-05
[21,] 0.99997465 5.069143e-05 2.534572e-05
[22,] 0.99996212 7.576049e-05 3.788024e-05
[23,] 0.99995138 9.723509e-05 4.861755e-05
[24,] 0.99990557 1.888662e-04 9.443308e-05
[25,] 0.99981888 3.622308e-04 1.811154e-04
[26,] 0.99970761 5.847810e-04 2.923905e-04
[27,] 0.99957233 8.553471e-04 4.276736e-04
[28,] 0.99950019 9.996266e-04 4.998133e-04
[29,] 0.99935368 1.292641e-03 6.463203e-04
[30,] 0.99897932 2.041358e-03 1.020679e-03
[31,] 0.99820386 3.592277e-03 1.796139e-03
[32,] 0.99700783 5.984338e-03 2.992169e-03
[33,] 0.99642305 7.153907e-03 3.576954e-03
[34,] 0.99713936 5.721287e-03 2.860644e-03
[35,] 0.99784251 4.314973e-03 2.157486e-03
[36,] 0.99778578 4.428439e-03 2.214219e-03
[37,] 0.99747719 5.045611e-03 2.522805e-03
[38,] 0.99723282 5.534357e-03 2.767179e-03
[39,] 0.99646621 7.067581e-03 3.533791e-03
[40,] 0.99500138 9.997240e-03 4.998620e-03
[41,] 0.99507681 9.846388e-03 4.923194e-03
[42,] 0.99550078 8.998434e-03 4.499217e-03
[43,] 0.99434662 1.130675e-02 5.653376e-03
[44,] 0.99336165 1.327669e-02 6.638346e-03
[45,] 0.99336716 1.326568e-02 6.632841e-03
[46,] 0.99551063 8.978734e-03 4.489367e-03
[47,] 0.99656321 6.873586e-03 3.436793e-03
[48,] 0.99769772 4.604554e-03 2.302277e-03
[49,] 0.99782405 4.351909e-03 2.175955e-03
[50,] 0.99704900 5.901994e-03 2.950997e-03
[51,] 0.99583027 8.339461e-03 4.169730e-03
[52,] 0.99363527 1.272946e-02 6.364729e-03
[53,] 0.99052590 1.894820e-02 9.474102e-03
[54,] 0.99576162 8.476763e-03 4.238381e-03
[55,] 0.99523463 9.530743e-03 4.765372e-03
[56,] 0.99349235 1.301530e-02 6.507650e-03
[57,] 0.99045199 1.909601e-02 9.548007e-03
[58,] 0.98568645 2.862709e-02 1.431355e-02
[59,] 0.97600184 4.799631e-02 2.399816e-02
[60,] 0.96934537 6.130927e-02 3.065463e-02
[61,] 0.95161516 9.676968e-02 4.838484e-02
[62,] 0.95325092 9.349816e-02 4.674908e-02
[63,] 0.98033970 3.932061e-02 1.966030e-02
> postscript(file="/var/www/html/rcomp/tmp/1e3kz1228386604.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/2kox01228386604.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/3bbu01228386605.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/4853z1228386605.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/5arzq1228386605.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 = 96
Frequency = 1
1 2 3 4 5 6
5.74310606 2.34560606 2.82560606 3.62060606 3.39685606 4.90310606
7 8 9 10 11 12
7.98185606 10.14310606 9.78560606 9.62176948 10.61801948 9.93676948
13 14 15 16 17 18
7.72793290 6.92043290 3.83043290 3.65543290 3.93168290 4.18793290
19 20 21 22 23 24
4.96668290 5.16793290 2.38043290 2.12659632 3.02284632 2.63159632
25 26 27 28 29 30
1.09275974 0.54525974 1.28525974 0.51025974 0.23650974 -2.25724026
31 32 33 34 35 36
-6.37849026 -7.46724026 -7.43474026 -7.97857684 -6.42232684 -8.02357684
37 38 39 40 41 42
-10.19241342 -11.44991342 -12.17991342 -10.47491342 -10.05866342 -8.91241342
43 44 45 46 47 48
-8.67366342 -8.42241342 -8.41991342 -9.46375000 -9.18750000 -8.69875000
49 50 51 52 53 54
-8.29758658 -6.94508658 -6.43508658 -6.93008658 -7.02383658 -5.90758658
55 56 57 58 59 60
-6.42883658 -6.71758658 -4.73508658 -4.73892316 -4.36267316 -4.13392316
61 62 63 64 65 66
-4.76275974 -4.05025974 -2.50025974 -3.48525974 -2.37900974 -1.17275974
67 68 69 70 71 72
-0.59400974 -0.69275974 -0.06025974 -0.63409632 0.29215368 1.95090368
73 74 75 76 77 78
1.95206710 4.29456710 4.96456710 3.96956710 4.11581710 2.23206710
79 80 81 82 83 84
2.97081710 5.14206710 7.30456710 8.43073052 7.56698052 7.44573052
85 86 87 88 89 90
6.73689394 8.33939394 8.20939394 9.13439394 7.78064394 6.92689394
91 92 93 94 95 96
6.15564394 2.84689394 1.17939394 2.63625000 -1.52750000 -1.10875000
> postscript(file="/var/www/html/rcomp/tmp/6qwdq1228386605.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 = 96
Frequency = 1
lag(myerror, k = 1) myerror
0 5.74310606 NA
1 2.34560606 5.74310606
2 2.82560606 2.34560606
3 3.62060606 2.82560606
4 3.39685606 3.62060606
5 4.90310606 3.39685606
6 7.98185606 4.90310606
7 10.14310606 7.98185606
8 9.78560606 10.14310606
9 9.62176948 9.78560606
10 10.61801948 9.62176948
11 9.93676948 10.61801948
12 7.72793290 9.93676948
13 6.92043290 7.72793290
14 3.83043290 6.92043290
15 3.65543290 3.83043290
16 3.93168290 3.65543290
17 4.18793290 3.93168290
18 4.96668290 4.18793290
19 5.16793290 4.96668290
20 2.38043290 5.16793290
21 2.12659632 2.38043290
22 3.02284632 2.12659632
23 2.63159632 3.02284632
24 1.09275974 2.63159632
25 0.54525974 1.09275974
26 1.28525974 0.54525974
27 0.51025974 1.28525974
28 0.23650974 0.51025974
29 -2.25724026 0.23650974
30 -6.37849026 -2.25724026
31 -7.46724026 -6.37849026
32 -7.43474026 -7.46724026
33 -7.97857684 -7.43474026
34 -6.42232684 -7.97857684
35 -8.02357684 -6.42232684
36 -10.19241342 -8.02357684
37 -11.44991342 -10.19241342
38 -12.17991342 -11.44991342
39 -10.47491342 -12.17991342
40 -10.05866342 -10.47491342
41 -8.91241342 -10.05866342
42 -8.67366342 -8.91241342
43 -8.42241342 -8.67366342
44 -8.41991342 -8.42241342
45 -9.46375000 -8.41991342
46 -9.18750000 -9.46375000
47 -8.69875000 -9.18750000
48 -8.29758658 -8.69875000
49 -6.94508658 -8.29758658
50 -6.43508658 -6.94508658
51 -6.93008658 -6.43508658
52 -7.02383658 -6.93008658
53 -5.90758658 -7.02383658
54 -6.42883658 -5.90758658
55 -6.71758658 -6.42883658
56 -4.73508658 -6.71758658
57 -4.73892316 -4.73508658
58 -4.36267316 -4.73892316
59 -4.13392316 -4.36267316
60 -4.76275974 -4.13392316
61 -4.05025974 -4.76275974
62 -2.50025974 -4.05025974
63 -3.48525974 -2.50025974
64 -2.37900974 -3.48525974
65 -1.17275974 -2.37900974
66 -0.59400974 -1.17275974
67 -0.69275974 -0.59400974
68 -0.06025974 -0.69275974
69 -0.63409632 -0.06025974
70 0.29215368 -0.63409632
71 1.95090368 0.29215368
72 1.95206710 1.95090368
73 4.29456710 1.95206710
74 4.96456710 4.29456710
75 3.96956710 4.96456710
76 4.11581710 3.96956710
77 2.23206710 4.11581710
78 2.97081710 2.23206710
79 5.14206710 2.97081710
80 7.30456710 5.14206710
81 8.43073052 7.30456710
82 7.56698052 8.43073052
83 7.44573052 7.56698052
84 6.73689394 7.44573052
85 8.33939394 6.73689394
86 8.20939394 8.33939394
87 9.13439394 8.20939394
88 7.78064394 9.13439394
89 6.92689394 7.78064394
90 6.15564394 6.92689394
91 2.84689394 6.15564394
92 1.17939394 2.84689394
93 2.63625000 1.17939394
94 -1.52750000 2.63625000
95 -1.10875000 -1.52750000
96 NA -1.10875000
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 2.34560606 5.74310606
[2,] 2.82560606 2.34560606
[3,] 3.62060606 2.82560606
[4,] 3.39685606 3.62060606
[5,] 4.90310606 3.39685606
[6,] 7.98185606 4.90310606
[7,] 10.14310606 7.98185606
[8,] 9.78560606 10.14310606
[9,] 9.62176948 9.78560606
[10,] 10.61801948 9.62176948
[11,] 9.93676948 10.61801948
[12,] 7.72793290 9.93676948
[13,] 6.92043290 7.72793290
[14,] 3.83043290 6.92043290
[15,] 3.65543290 3.83043290
[16,] 3.93168290 3.65543290
[17,] 4.18793290 3.93168290
[18,] 4.96668290 4.18793290
[19,] 5.16793290 4.96668290
[20,] 2.38043290 5.16793290
[21,] 2.12659632 2.38043290
[22,] 3.02284632 2.12659632
[23,] 2.63159632 3.02284632
[24,] 1.09275974 2.63159632
[25,] 0.54525974 1.09275974
[26,] 1.28525974 0.54525974
[27,] 0.51025974 1.28525974
[28,] 0.23650974 0.51025974
[29,] -2.25724026 0.23650974
[30,] -6.37849026 -2.25724026
[31,] -7.46724026 -6.37849026
[32,] -7.43474026 -7.46724026
[33,] -7.97857684 -7.43474026
[34,] -6.42232684 -7.97857684
[35,] -8.02357684 -6.42232684
[36,] -10.19241342 -8.02357684
[37,] -11.44991342 -10.19241342
[38,] -12.17991342 -11.44991342
[39,] -10.47491342 -12.17991342
[40,] -10.05866342 -10.47491342
[41,] -8.91241342 -10.05866342
[42,] -8.67366342 -8.91241342
[43,] -8.42241342 -8.67366342
[44,] -8.41991342 -8.42241342
[45,] -9.46375000 -8.41991342
[46,] -9.18750000 -9.46375000
[47,] -8.69875000 -9.18750000
[48,] -8.29758658 -8.69875000
[49,] -6.94508658 -8.29758658
[50,] -6.43508658 -6.94508658
[51,] -6.93008658 -6.43508658
[52,] -7.02383658 -6.93008658
[53,] -5.90758658 -7.02383658
[54,] -6.42883658 -5.90758658
[55,] -6.71758658 -6.42883658
[56,] -4.73508658 -6.71758658
[57,] -4.73892316 -4.73508658
[58,] -4.36267316 -4.73892316
[59,] -4.13392316 -4.36267316
[60,] -4.76275974 -4.13392316
[61,] -4.05025974 -4.76275974
[62,] -2.50025974 -4.05025974
[63,] -3.48525974 -2.50025974
[64,] -2.37900974 -3.48525974
[65,] -1.17275974 -2.37900974
[66,] -0.59400974 -1.17275974
[67,] -0.69275974 -0.59400974
[68,] -0.06025974 -0.69275974
[69,] -0.63409632 -0.06025974
[70,] 0.29215368 -0.63409632
[71,] 1.95090368 0.29215368
[72,] 1.95206710 1.95090368
[73,] 4.29456710 1.95206710
[74,] 4.96456710 4.29456710
[75,] 3.96956710 4.96456710
[76,] 4.11581710 3.96956710
[77,] 2.23206710 4.11581710
[78,] 2.97081710 2.23206710
[79,] 5.14206710 2.97081710
[80,] 7.30456710 5.14206710
[81,] 8.43073052 7.30456710
[82,] 7.56698052 8.43073052
[83,] 7.44573052 7.56698052
[84,] 6.73689394 7.44573052
[85,] 8.33939394 6.73689394
[86,] 8.20939394 8.33939394
[87,] 9.13439394 8.20939394
[88,] 7.78064394 9.13439394
[89,] 6.92689394 7.78064394
[90,] 6.15564394 6.92689394
[91,] 2.84689394 6.15564394
[92,] 1.17939394 2.84689394
[93,] 2.63625000 1.17939394
[94,] -1.52750000 2.63625000
[95,] -1.10875000 -1.52750000
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 2.34560606 5.74310606
2 2.82560606 2.34560606
3 3.62060606 2.82560606
4 3.39685606 3.62060606
5 4.90310606 3.39685606
6 7.98185606 4.90310606
7 10.14310606 7.98185606
8 9.78560606 10.14310606
9 9.62176948 9.78560606
10 10.61801948 9.62176948
11 9.93676948 10.61801948
12 7.72793290 9.93676948
13 6.92043290 7.72793290
14 3.83043290 6.92043290
15 3.65543290 3.83043290
16 3.93168290 3.65543290
17 4.18793290 3.93168290
18 4.96668290 4.18793290
19 5.16793290 4.96668290
20 2.38043290 5.16793290
21 2.12659632 2.38043290
22 3.02284632 2.12659632
23 2.63159632 3.02284632
24 1.09275974 2.63159632
25 0.54525974 1.09275974
26 1.28525974 0.54525974
27 0.51025974 1.28525974
28 0.23650974 0.51025974
29 -2.25724026 0.23650974
30 -6.37849026 -2.25724026
31 -7.46724026 -6.37849026
32 -7.43474026 -7.46724026
33 -7.97857684 -7.43474026
34 -6.42232684 -7.97857684
35 -8.02357684 -6.42232684
36 -10.19241342 -8.02357684
37 -11.44991342 -10.19241342
38 -12.17991342 -11.44991342
39 -10.47491342 -12.17991342
40 -10.05866342 -10.47491342
41 -8.91241342 -10.05866342
42 -8.67366342 -8.91241342
43 -8.42241342 -8.67366342
44 -8.41991342 -8.42241342
45 -9.46375000 -8.41991342
46 -9.18750000 -9.46375000
47 -8.69875000 -9.18750000
48 -8.29758658 -8.69875000
49 -6.94508658 -8.29758658
50 -6.43508658 -6.94508658
51 -6.93008658 -6.43508658
52 -7.02383658 -6.93008658
53 -5.90758658 -7.02383658
54 -6.42883658 -5.90758658
55 -6.71758658 -6.42883658
56 -4.73508658 -6.71758658
57 -4.73892316 -4.73508658
58 -4.36267316 -4.73892316
59 -4.13392316 -4.36267316
60 -4.76275974 -4.13392316
61 -4.05025974 -4.76275974
62 -2.50025974 -4.05025974
63 -3.48525974 -2.50025974
64 -2.37900974 -3.48525974
65 -1.17275974 -2.37900974
66 -0.59400974 -1.17275974
67 -0.69275974 -0.59400974
68 -0.06025974 -0.69275974
69 -0.63409632 -0.06025974
70 0.29215368 -0.63409632
71 1.95090368 0.29215368
72 1.95206710 1.95090368
73 4.29456710 1.95206710
74 4.96456710 4.29456710
75 3.96956710 4.96456710
76 4.11581710 3.96956710
77 2.23206710 4.11581710
78 2.97081710 2.23206710
79 5.14206710 2.97081710
80 7.30456710 5.14206710
81 8.43073052 7.30456710
82 7.56698052 8.43073052
83 7.44573052 7.56698052
84 6.73689394 7.44573052
85 8.33939394 6.73689394
86 8.20939394 8.33939394
87 9.13439394 8.20939394
88 7.78064394 9.13439394
89 6.92689394 7.78064394
90 6.15564394 6.92689394
91 2.84689394 6.15564394
92 1.17939394 2.84689394
93 2.63625000 1.17939394
94 -1.52750000 2.63625000
95 -1.10875000 -1.52750000
> 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/7rgtl1228386605.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/81s601228386605.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/96a381228386605.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/10i16c1228386605.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/11bzyb1228386605.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/12e9it1228386605.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/13zqva1228386605.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/14aa0x1228386605.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/152fea1228386605.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/161sxp1228386605.tab")
+ }
>
> system("convert tmp/1e3kz1228386604.ps tmp/1e3kz1228386604.png")
> system("convert tmp/2kox01228386604.ps tmp/2kox01228386604.png")
> system("convert tmp/3bbu01228386605.ps tmp/3bbu01228386605.png")
> system("convert tmp/4853z1228386605.ps tmp/4853z1228386605.png")
> system("convert tmp/5arzq1228386605.ps tmp/5arzq1228386605.png")
> system("convert tmp/6qwdq1228386605.ps tmp/6qwdq1228386605.png")
> system("convert tmp/7rgtl1228386605.ps tmp/7rgtl1228386605.png")
> system("convert tmp/81s601228386605.ps tmp/81s601228386605.png")
> system("convert tmp/96a381228386605.ps tmp/96a381228386605.png")
> system("convert tmp/10i16c1228386605.ps tmp/10i16c1228386605.png")
>
>
> proc.time()
user system elapsed
2.937 1.645 3.881