R version 2.7.0 (2008-04-22)
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(5.5,0,5.3,0,5.2,0,5.3,0,5.3,0,5,0,4.8,0,4.9,0,5.3,0,6,0,6.2,0,6.4,0,6.4,0,6.4,0,6.2,0,6.1,0,6,0,5.9,0,6.2,0,6.2,0,6.4,0,6.8,0,6.9,0,7,0,7,1,6.9,1,6.7,1,6.6,1,6.5,1,6.4,1,6.5,1,6.5,1,6.6,1,6.7,1,6.8,1,7.2,1,7.6,1,7.6,1,7.3,1,6.4,1,6.1,1,6.3,1,7.1,1,7.5,1,7.4,1,7.1,1,6.8,1,6.9,1,7.2,1,7.4,1,7.3,1,6.9,1,6.9,1,6.8,1,7.1,1,7.2,1,7.1,1,7,1,6.9,1,7,1,7.4,1,7.5,1,7.5,1,7.4,1,7.3,1,7,1,6.7,1,6.5,1,6.5,1,6.5,1,6.6,1,6.8,1,6.9,1,6.9,1,6.8,1,6.8,1,6.5,1,6.1,1,6,1,5.9,1,5.8,1,5.9,1,5.9,1,6.2,1,6.3,1,6.2,1,6,1,5.8,1,5.5,1,5.5,1,5.7,1,5.8,1),dim=c(2,92),dimnames=list(c('VAR1','D1'),1:92))
> y <- array(NA,dim=c(2,92),dimnames=list(c('VAR1','D1'),1:92))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par20 = ''
> par19 = ''
> par18 = ''
> par17 = ''
> par16 = ''
> par15 = ''
> par14 = ''
> par13 = ''
> par12 = ''
> par11 = ''
> par10 = ''
> par9 = ''
> par8 = ''
> par7 = ''
> par6 = ''
> par5 = ''
> par4 = ''
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> ylab = ''
> xlab = ''
> main = ''
> #'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
VAR1 D1 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 5.5 0 1 0 0 0 0 0 0 0 0 0 0 1
2 5.3 0 0 1 0 0 0 0 0 0 0 0 0 2
3 5.2 0 0 0 1 0 0 0 0 0 0 0 0 3
4 5.3 0 0 0 0 1 0 0 0 0 0 0 0 4
5 5.3 0 0 0 0 0 1 0 0 0 0 0 0 5
6 5.0 0 0 0 0 0 0 1 0 0 0 0 0 6
7 4.8 0 0 0 0 0 0 0 1 0 0 0 0 7
8 4.9 0 0 0 0 0 0 0 0 1 0 0 0 8
9 5.3 0 0 0 0 0 0 0 0 0 1 0 0 9
10 6.0 0 0 0 0 0 0 0 0 0 0 1 0 10
11 6.2 0 0 0 0 0 0 0 0 0 0 0 1 11
12 6.4 0 0 0 0 0 0 0 0 0 0 0 0 12
13 6.4 0 1 0 0 0 0 0 0 0 0 0 0 13
14 6.4 0 0 1 0 0 0 0 0 0 0 0 0 14
15 6.2 0 0 0 1 0 0 0 0 0 0 0 0 15
16 6.1 0 0 0 0 1 0 0 0 0 0 0 0 16
17 6.0 0 0 0 0 0 1 0 0 0 0 0 0 17
18 5.9 0 0 0 0 0 0 1 0 0 0 0 0 18
19 6.2 0 0 0 0 0 0 0 1 0 0 0 0 19
20 6.2 0 0 0 0 0 0 0 0 1 0 0 0 20
21 6.4 0 0 0 0 0 0 0 0 0 1 0 0 21
22 6.8 0 0 0 0 0 0 0 0 0 0 1 0 22
23 6.9 0 0 0 0 0 0 0 0 0 0 0 1 23
24 7.0 0 0 0 0 0 0 0 0 0 0 0 0 24
25 7.0 1 1 0 0 0 0 0 0 0 0 0 0 25
26 6.9 1 0 1 0 0 0 0 0 0 0 0 0 26
27 6.7 1 0 0 1 0 0 0 0 0 0 0 0 27
28 6.6 1 0 0 0 1 0 0 0 0 0 0 0 28
29 6.5 1 0 0 0 0 1 0 0 0 0 0 0 29
30 6.4 1 0 0 0 0 0 1 0 0 0 0 0 30
31 6.5 1 0 0 0 0 0 0 1 0 0 0 0 31
32 6.5 1 0 0 0 0 0 0 0 1 0 0 0 32
33 6.6 1 0 0 0 0 0 0 0 0 1 0 0 33
34 6.7 1 0 0 0 0 0 0 0 0 0 1 0 34
35 6.8 1 0 0 0 0 0 0 0 0 0 0 1 35
36 7.2 1 0 0 0 0 0 0 0 0 0 0 0 36
37 7.6 1 1 0 0 0 0 0 0 0 0 0 0 37
38 7.6 1 0 1 0 0 0 0 0 0 0 0 0 38
39 7.3 1 0 0 1 0 0 0 0 0 0 0 0 39
40 6.4 1 0 0 0 1 0 0 0 0 0 0 0 40
41 6.1 1 0 0 0 0 1 0 0 0 0 0 0 41
42 6.3 1 0 0 0 0 0 1 0 0 0 0 0 42
43 7.1 1 0 0 0 0 0 0 1 0 0 0 0 43
44 7.5 1 0 0 0 0 0 0 0 1 0 0 0 44
45 7.4 1 0 0 0 0 0 0 0 0 1 0 0 45
46 7.1 1 0 0 0 0 0 0 0 0 0 1 0 46
47 6.8 1 0 0 0 0 0 0 0 0 0 0 1 47
48 6.9 1 0 0 0 0 0 0 0 0 0 0 0 48
49 7.2 1 1 0 0 0 0 0 0 0 0 0 0 49
50 7.4 1 0 1 0 0 0 0 0 0 0 0 0 50
51 7.3 1 0 0 1 0 0 0 0 0 0 0 0 51
52 6.9 1 0 0 0 1 0 0 0 0 0 0 0 52
53 6.9 1 0 0 0 0 1 0 0 0 0 0 0 53
54 6.8 1 0 0 0 0 0 1 0 0 0 0 0 54
55 7.1 1 0 0 0 0 0 0 1 0 0 0 0 55
56 7.2 1 0 0 0 0 0 0 0 1 0 0 0 56
57 7.1 1 0 0 0 0 0 0 0 0 1 0 0 57
58 7.0 1 0 0 0 0 0 0 0 0 0 1 0 58
59 6.9 1 0 0 0 0 0 0 0 0 0 0 1 59
60 7.0 1 0 0 0 0 0 0 0 0 0 0 0 60
61 7.4 1 1 0 0 0 0 0 0 0 0 0 0 61
62 7.5 1 0 1 0 0 0 0 0 0 0 0 0 62
63 7.5 1 0 0 1 0 0 0 0 0 0 0 0 63
64 7.4 1 0 0 0 1 0 0 0 0 0 0 0 64
65 7.3 1 0 0 0 0 1 0 0 0 0 0 0 65
66 7.0 1 0 0 0 0 0 1 0 0 0 0 0 66
67 6.7 1 0 0 0 0 0 0 1 0 0 0 0 67
68 6.5 1 0 0 0 0 0 0 0 1 0 0 0 68
69 6.5 1 0 0 0 0 0 0 0 0 1 0 0 69
70 6.5 1 0 0 0 0 0 0 0 0 0 1 0 70
71 6.6 1 0 0 0 0 0 0 0 0 0 0 1 71
72 6.8 1 0 0 0 0 0 0 0 0 0 0 0 72
73 6.9 1 1 0 0 0 0 0 0 0 0 0 0 73
74 6.9 1 0 1 0 0 0 0 0 0 0 0 0 74
75 6.8 1 0 0 1 0 0 0 0 0 0 0 0 75
76 6.8 1 0 0 0 1 0 0 0 0 0 0 0 76
77 6.5 1 0 0 0 0 1 0 0 0 0 0 0 77
78 6.1 1 0 0 0 0 0 1 0 0 0 0 0 78
79 6.0 1 0 0 0 0 0 0 1 0 0 0 0 79
80 5.9 1 0 0 0 0 0 0 0 1 0 0 0 80
81 5.8 1 0 0 0 0 0 0 0 0 1 0 0 81
82 5.9 1 0 0 0 0 0 0 0 0 0 1 0 82
83 5.9 1 0 0 0 0 0 0 0 0 0 0 1 83
84 6.2 1 0 0 0 0 0 0 0 0 0 0 0 84
85 6.3 1 1 0 0 0 0 0 0 0 0 0 0 85
86 6.2 1 0 1 0 0 0 0 0 0 0 0 0 86
87 6.0 1 0 0 1 0 0 0 0 0 0 0 0 87
88 5.8 1 0 0 0 1 0 0 0 0 0 0 0 88
89 5.5 1 0 0 0 0 1 0 0 0 0 0 0 89
90 5.5 1 0 0 0 0 0 1 0 0 0 0 0 90
91 5.7 1 0 0 0 0 0 0 1 0 0 0 0 91
92 5.8 1 0 0 0 0 0 0 0 1 0 0 0 92
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) D1 M1 M2 M3 M4
6.36969 1.27868 -0.09569 -0.09783 -0.23746 -0.43960
M5 M6 M7 M8 M9 M10
-0.57924 -0.70638 -0.55852 -0.49816 -0.37394 -0.23501
M11 t
-0.21036 -0.01036
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.9511 -0.3950 -0.0289 0.3933 0.9790
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.369695 0.225141 28.292 < 2e-16 ***
D1 1.278684 0.195103 6.554 5.49e-09 ***
M1 -0.095686 0.275426 -0.347 0.72922
M2 -0.097825 0.275192 -0.355 0.72319
M3 -0.237464 0.274997 -0.864 0.39050
M4 -0.439603 0.274839 -1.599 0.11375
M5 -0.579243 0.274720 -2.108 0.03820 *
M6 -0.706382 0.274638 -2.572 0.01201 *
M7 -0.558521 0.274595 -2.034 0.04535 *
M8 -0.498160 0.274589 -1.814 0.07349 .
M9 -0.373940 0.283720 -1.318 0.19136
M10 -0.235008 0.283627 -0.829 0.40987
M11 -0.210361 0.283572 -0.742 0.46042
t -0.010361 0.003233 -3.205 0.00196 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5305 on 78 degrees of freedom
Multiple R-squared: 0.4601, Adjusted R-squared: 0.3701
F-statistic: 5.112 on 13 and 78 DF, p-value: 1.826e-06
> 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,] 3.233281e-02 6.466561e-02 0.9676671938
[2,] 7.718303e-03 1.543661e-02 0.9922816966
[3,] 2.373105e-02 4.746210e-02 0.9762689517
[4,] 1.688853e-02 3.377706e-02 0.9831114695
[5,] 6.674011e-03 1.334802e-02 0.9933259887
[6,] 3.266419e-03 6.532837e-03 0.9967335813
[7,] 2.107063e-03 4.214126e-03 0.9978929370
[8,] 1.811059e-03 3.622118e-03 0.9981889411
[9,] 7.464898e-04 1.492980e-03 0.9992535102
[10,] 3.314743e-04 6.629487e-04 0.9996685257
[11,] 1.629840e-04 3.259680e-04 0.9998370160
[12,] 8.215326e-05 1.643065e-04 0.9999178467
[13,] 4.266518e-05 8.533035e-05 0.9999573348
[14,] 1.828944e-05 3.657889e-05 0.9999817106
[15,] 8.945831e-06 1.789166e-05 0.9999910542
[16,] 5.185848e-06 1.037170e-05 0.9999948142
[17,] 5.272778e-06 1.054556e-05 0.9999947272
[18,] 1.443095e-04 2.886190e-04 0.9998556905
[19,] 8.800972e-04 1.760194e-03 0.9991199028
[20,] 7.723758e-04 1.544752e-03 0.9992276242
[21,] 4.428534e-04 8.857067e-04 0.9995571466
[22,] 2.185351e-04 4.370702e-04 0.9997814649
[23,] 1.416929e-04 2.833857e-04 0.9998583071
[24,] 1.801338e-02 3.602675e-02 0.9819866230
[25,] 3.841346e-01 7.682691e-01 0.6158654347
[26,] 7.036875e-01 5.926250e-01 0.2963125125
[27,] 6.687385e-01 6.625230e-01 0.3312615168
[28,] 6.577950e-01 6.844100e-01 0.3422050219
[29,] 5.950309e-01 8.099381e-01 0.4049690711
[30,] 6.243489e-01 7.513023e-01 0.3756511370
[31,] 8.093511e-01 3.812979e-01 0.1906489478
[32,] 9.389633e-01 1.220735e-01 0.0610367489
[33,] 9.731032e-01 5.379350e-02 0.0268967502
[34,] 9.764926e-01 4.701489e-02 0.0235074468
[35,] 9.794484e-01 4.110328e-02 0.0205516393
[36,] 9.963096e-01 7.380773e-03 0.0036903867
[37,] 9.990545e-01 1.891047e-03 0.0009455234
[38,] 9.997617e-01 4.766432e-04 0.0002383216
[39,] 9.996757e-01 6.486110e-04 0.0003243055
[40,] 9.993788e-01 1.242408e-03 0.0006212041
[41,] 9.989315e-01 2.137069e-03 0.0010685346
[42,] 9.987263e-01 2.547343e-03 0.0012736713
[43,] 9.990258e-01 1.948484e-03 0.0009742421
[44,] 9.997212e-01 5.576326e-04 0.0002788163
[45,] 9.995670e-01 8.660006e-04 0.0004330003
[46,] 9.990679e-01 1.864186e-03 0.0009320931
[47,] 9.980791e-01 3.841756e-03 0.0019208781
[48,] 9.962455e-01 7.509042e-03 0.0037545212
[49,] 9.963488e-01 7.302380e-03 0.0036511898
[50,] 9.956090e-01 8.782016e-03 0.0043910078
[51,] 9.921837e-01 1.563256e-02 0.0078162791
[52,] 9.929350e-01 1.412998e-02 0.0070649882
[53,] 9.885054e-01 2.298917e-02 0.0114945837
[54,] 9.812924e-01 3.741530e-02 0.0187076478
[55,] 9.663814e-01 6.723720e-02 0.0336186010
[56,] 9.375292e-01 1.249417e-01 0.0624708428
[57,] 8.829172e-01 2.341657e-01 0.1170828457
[58,] 7.898949e-01 4.202103e-01 0.2101051271
[59,] 6.581625e-01 6.836751e-01 0.3418375281
> postscript(file="/var/www/html/rcomp/tmp/1hsxq1229610378.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/28h191229610378.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/3sqj11229610378.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/49tzv1229610378.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/5ck241229610378.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 = 92
Frequency = 1
1 2 3 4 5 6
-0.76364734 -0.95114734 -0.90114734 -0.58864734 -0.43864734 -0.60114734
7 8 9 10 11 12
-0.93864734 -0.88864734 -0.60250604 -0.03107747 0.15463682 0.15463682
13 14 15 16 17 18
0.26068409 0.27318409 0.22318409 0.33568409 0.38568409 0.42318409
19 20 21 22 23 24
0.58568409 0.53568409 0.62182540 0.89325397 0.97896825 0.87896825
25 26 27 28 29 30
-0.29366805 -0.38116805 -0.43116805 -0.31866805 -0.26866805 -0.23116805
31 32 33 34 35 36
-0.26866805 -0.31866805 -0.33252674 -0.36109817 -0.27538389 -0.07538389
37 38 39 40 41 42
0.43066339 0.44316339 0.29316339 -0.39433661 -0.54433661 -0.20683661
43 44 45 46 47 48
0.45566339 0.80566339 0.59180469 0.16323326 -0.15105245 -0.25105245
49 50 51 52 53 54
0.15499482 0.36749482 0.41749482 0.22999482 0.37999482 0.41749482
55 56 57 58 59 60
0.57999482 0.62999482 0.41613613 0.18756470 0.07327899 -0.02672101
61 62 63 64 65 66
0.47932626 0.59182626 0.74182626 0.85432626 0.90432626 0.74182626
67 68 69 70 71 72
0.30432626 0.05432626 -0.05953244 -0.18810386 -0.10238958 -0.10238958
73 74 75 76 77 78
0.10365769 0.11615769 0.16615769 0.37865769 0.22865769 -0.03384231
79 80 81 82 83 84
-0.27134231 -0.42134231 -0.63520100 -0.66377243 -0.67805814 -0.57805814
85 86 87 88 89 90
-0.37201087 -0.45951087 -0.50951087 -0.49701087 -0.64701087 -0.50951087
91 92
-0.44701087 -0.39701087
> postscript(file="/var/www/html/rcomp/tmp/6sjk41229610378.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 = 92
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.76364734 NA
1 -0.95114734 -0.76364734
2 -0.90114734 -0.95114734
3 -0.58864734 -0.90114734
4 -0.43864734 -0.58864734
5 -0.60114734 -0.43864734
6 -0.93864734 -0.60114734
7 -0.88864734 -0.93864734
8 -0.60250604 -0.88864734
9 -0.03107747 -0.60250604
10 0.15463682 -0.03107747
11 0.15463682 0.15463682
12 0.26068409 0.15463682
13 0.27318409 0.26068409
14 0.22318409 0.27318409
15 0.33568409 0.22318409
16 0.38568409 0.33568409
17 0.42318409 0.38568409
18 0.58568409 0.42318409
19 0.53568409 0.58568409
20 0.62182540 0.53568409
21 0.89325397 0.62182540
22 0.97896825 0.89325397
23 0.87896825 0.97896825
24 -0.29366805 0.87896825
25 -0.38116805 -0.29366805
26 -0.43116805 -0.38116805
27 -0.31866805 -0.43116805
28 -0.26866805 -0.31866805
29 -0.23116805 -0.26866805
30 -0.26866805 -0.23116805
31 -0.31866805 -0.26866805
32 -0.33252674 -0.31866805
33 -0.36109817 -0.33252674
34 -0.27538389 -0.36109817
35 -0.07538389 -0.27538389
36 0.43066339 -0.07538389
37 0.44316339 0.43066339
38 0.29316339 0.44316339
39 -0.39433661 0.29316339
40 -0.54433661 -0.39433661
41 -0.20683661 -0.54433661
42 0.45566339 -0.20683661
43 0.80566339 0.45566339
44 0.59180469 0.80566339
45 0.16323326 0.59180469
46 -0.15105245 0.16323326
47 -0.25105245 -0.15105245
48 0.15499482 -0.25105245
49 0.36749482 0.15499482
50 0.41749482 0.36749482
51 0.22999482 0.41749482
52 0.37999482 0.22999482
53 0.41749482 0.37999482
54 0.57999482 0.41749482
55 0.62999482 0.57999482
56 0.41613613 0.62999482
57 0.18756470 0.41613613
58 0.07327899 0.18756470
59 -0.02672101 0.07327899
60 0.47932626 -0.02672101
61 0.59182626 0.47932626
62 0.74182626 0.59182626
63 0.85432626 0.74182626
64 0.90432626 0.85432626
65 0.74182626 0.90432626
66 0.30432626 0.74182626
67 0.05432626 0.30432626
68 -0.05953244 0.05432626
69 -0.18810386 -0.05953244
70 -0.10238958 -0.18810386
71 -0.10238958 -0.10238958
72 0.10365769 -0.10238958
73 0.11615769 0.10365769
74 0.16615769 0.11615769
75 0.37865769 0.16615769
76 0.22865769 0.37865769
77 -0.03384231 0.22865769
78 -0.27134231 -0.03384231
79 -0.42134231 -0.27134231
80 -0.63520100 -0.42134231
81 -0.66377243 -0.63520100
82 -0.67805814 -0.66377243
83 -0.57805814 -0.67805814
84 -0.37201087 -0.57805814
85 -0.45951087 -0.37201087
86 -0.50951087 -0.45951087
87 -0.49701087 -0.50951087
88 -0.64701087 -0.49701087
89 -0.50951087 -0.64701087
90 -0.44701087 -0.50951087
91 -0.39701087 -0.44701087
92 NA -0.39701087
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.95114734 -0.76364734
[2,] -0.90114734 -0.95114734
[3,] -0.58864734 -0.90114734
[4,] -0.43864734 -0.58864734
[5,] -0.60114734 -0.43864734
[6,] -0.93864734 -0.60114734
[7,] -0.88864734 -0.93864734
[8,] -0.60250604 -0.88864734
[9,] -0.03107747 -0.60250604
[10,] 0.15463682 -0.03107747
[11,] 0.15463682 0.15463682
[12,] 0.26068409 0.15463682
[13,] 0.27318409 0.26068409
[14,] 0.22318409 0.27318409
[15,] 0.33568409 0.22318409
[16,] 0.38568409 0.33568409
[17,] 0.42318409 0.38568409
[18,] 0.58568409 0.42318409
[19,] 0.53568409 0.58568409
[20,] 0.62182540 0.53568409
[21,] 0.89325397 0.62182540
[22,] 0.97896825 0.89325397
[23,] 0.87896825 0.97896825
[24,] -0.29366805 0.87896825
[25,] -0.38116805 -0.29366805
[26,] -0.43116805 -0.38116805
[27,] -0.31866805 -0.43116805
[28,] -0.26866805 -0.31866805
[29,] -0.23116805 -0.26866805
[30,] -0.26866805 -0.23116805
[31,] -0.31866805 -0.26866805
[32,] -0.33252674 -0.31866805
[33,] -0.36109817 -0.33252674
[34,] -0.27538389 -0.36109817
[35,] -0.07538389 -0.27538389
[36,] 0.43066339 -0.07538389
[37,] 0.44316339 0.43066339
[38,] 0.29316339 0.44316339
[39,] -0.39433661 0.29316339
[40,] -0.54433661 -0.39433661
[41,] -0.20683661 -0.54433661
[42,] 0.45566339 -0.20683661
[43,] 0.80566339 0.45566339
[44,] 0.59180469 0.80566339
[45,] 0.16323326 0.59180469
[46,] -0.15105245 0.16323326
[47,] -0.25105245 -0.15105245
[48,] 0.15499482 -0.25105245
[49,] 0.36749482 0.15499482
[50,] 0.41749482 0.36749482
[51,] 0.22999482 0.41749482
[52,] 0.37999482 0.22999482
[53,] 0.41749482 0.37999482
[54,] 0.57999482 0.41749482
[55,] 0.62999482 0.57999482
[56,] 0.41613613 0.62999482
[57,] 0.18756470 0.41613613
[58,] 0.07327899 0.18756470
[59,] -0.02672101 0.07327899
[60,] 0.47932626 -0.02672101
[61,] 0.59182626 0.47932626
[62,] 0.74182626 0.59182626
[63,] 0.85432626 0.74182626
[64,] 0.90432626 0.85432626
[65,] 0.74182626 0.90432626
[66,] 0.30432626 0.74182626
[67,] 0.05432626 0.30432626
[68,] -0.05953244 0.05432626
[69,] -0.18810386 -0.05953244
[70,] -0.10238958 -0.18810386
[71,] -0.10238958 -0.10238958
[72,] 0.10365769 -0.10238958
[73,] 0.11615769 0.10365769
[74,] 0.16615769 0.11615769
[75,] 0.37865769 0.16615769
[76,] 0.22865769 0.37865769
[77,] -0.03384231 0.22865769
[78,] -0.27134231 -0.03384231
[79,] -0.42134231 -0.27134231
[80,] -0.63520100 -0.42134231
[81,] -0.66377243 -0.63520100
[82,] -0.67805814 -0.66377243
[83,] -0.57805814 -0.67805814
[84,] -0.37201087 -0.57805814
[85,] -0.45951087 -0.37201087
[86,] -0.50951087 -0.45951087
[87,] -0.49701087 -0.50951087
[88,] -0.64701087 -0.49701087
[89,] -0.50951087 -0.64701087
[90,] -0.44701087 -0.50951087
[91,] -0.39701087 -0.44701087
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.95114734 -0.76364734
2 -0.90114734 -0.95114734
3 -0.58864734 -0.90114734
4 -0.43864734 -0.58864734
5 -0.60114734 -0.43864734
6 -0.93864734 -0.60114734
7 -0.88864734 -0.93864734
8 -0.60250604 -0.88864734
9 -0.03107747 -0.60250604
10 0.15463682 -0.03107747
11 0.15463682 0.15463682
12 0.26068409 0.15463682
13 0.27318409 0.26068409
14 0.22318409 0.27318409
15 0.33568409 0.22318409
16 0.38568409 0.33568409
17 0.42318409 0.38568409
18 0.58568409 0.42318409
19 0.53568409 0.58568409
20 0.62182540 0.53568409
21 0.89325397 0.62182540
22 0.97896825 0.89325397
23 0.87896825 0.97896825
24 -0.29366805 0.87896825
25 -0.38116805 -0.29366805
26 -0.43116805 -0.38116805
27 -0.31866805 -0.43116805
28 -0.26866805 -0.31866805
29 -0.23116805 -0.26866805
30 -0.26866805 -0.23116805
31 -0.31866805 -0.26866805
32 -0.33252674 -0.31866805
33 -0.36109817 -0.33252674
34 -0.27538389 -0.36109817
35 -0.07538389 -0.27538389
36 0.43066339 -0.07538389
37 0.44316339 0.43066339
38 0.29316339 0.44316339
39 -0.39433661 0.29316339
40 -0.54433661 -0.39433661
41 -0.20683661 -0.54433661
42 0.45566339 -0.20683661
43 0.80566339 0.45566339
44 0.59180469 0.80566339
45 0.16323326 0.59180469
46 -0.15105245 0.16323326
47 -0.25105245 -0.15105245
48 0.15499482 -0.25105245
49 0.36749482 0.15499482
50 0.41749482 0.36749482
51 0.22999482 0.41749482
52 0.37999482 0.22999482
53 0.41749482 0.37999482
54 0.57999482 0.41749482
55 0.62999482 0.57999482
56 0.41613613 0.62999482
57 0.18756470 0.41613613
58 0.07327899 0.18756470
59 -0.02672101 0.07327899
60 0.47932626 -0.02672101
61 0.59182626 0.47932626
62 0.74182626 0.59182626
63 0.85432626 0.74182626
64 0.90432626 0.85432626
65 0.74182626 0.90432626
66 0.30432626 0.74182626
67 0.05432626 0.30432626
68 -0.05953244 0.05432626
69 -0.18810386 -0.05953244
70 -0.10238958 -0.18810386
71 -0.10238958 -0.10238958
72 0.10365769 -0.10238958
73 0.11615769 0.10365769
74 0.16615769 0.11615769
75 0.37865769 0.16615769
76 0.22865769 0.37865769
77 -0.03384231 0.22865769
78 -0.27134231 -0.03384231
79 -0.42134231 -0.27134231
80 -0.63520100 -0.42134231
81 -0.66377243 -0.63520100
82 -0.67805814 -0.66377243
83 -0.57805814 -0.67805814
84 -0.37201087 -0.57805814
85 -0.45951087 -0.37201087
86 -0.50951087 -0.45951087
87 -0.49701087 -0.50951087
88 -0.64701087 -0.49701087
89 -0.50951087 -0.64701087
90 -0.44701087 -0.50951087
91 -0.39701087 -0.44701087
> 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/7nmmy1229610378.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/86tq21229610378.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/9ej8b1229610378.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/10sxny1229610378.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/11e0fj1229610378.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/12481o1229610379.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/13cokw1229610379.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/14vg391229610379.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/15lsqt1229610379.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/161vw01229610379.tab")
+ }
>
> system("convert tmp/1hsxq1229610378.ps tmp/1hsxq1229610378.png")
> system("convert tmp/28h191229610378.ps tmp/28h191229610378.png")
> system("convert tmp/3sqj11229610378.ps tmp/3sqj11229610378.png")
> system("convert tmp/49tzv1229610378.ps tmp/49tzv1229610378.png")
> system("convert tmp/5ck241229610378.ps tmp/5ck241229610378.png")
> system("convert tmp/6sjk41229610378.ps tmp/6sjk41229610378.png")
> system("convert tmp/7nmmy1229610378.ps tmp/7nmmy1229610378.png")
> system("convert tmp/86tq21229610378.ps tmp/86tq21229610378.png")
> system("convert tmp/9ej8b1229610378.ps tmp/9ej8b1229610378.png")
> system("convert tmp/10sxny1229610378.ps tmp/10sxny1229610378.png")
>
>
> proc.time()
user system elapsed
5.678 2.806 6.059