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(96.8,0,87.0,0,96.3,0,107.1,0,115.2,0,106.1,0,89.5,0,91.3,0,97.6,0,100.7,0,104.6,0,94.7,0,101.8,0,102.5,0,105.3,0,110.3,0,109.8,0,117.3,0,118.8,0,131.3,0,125.9,0,133.1,0,147.0,0,145.8,0,164.4,0,149.8,0,137.7,0,151.7,0,156.8,0,180.0,0,180.4,0,170.4,0,191.6,0,199.5,0,218.2,0,217.5,0,205.0,0,194.0,0,199.3,0,219.3,0,211.1,0,215.2,0,240.2,0,242.2,0,240.7,0,255.4,0,253.0,0,218.2,0,203.7,0,205.6,0,215.6,0,188.5,0,202.9,0,214.0,0,230.3,0,230.0,0,241.0,0,259.6,1,247.8,1,270.3,1,289.7,1,322.7,1,315.0,1,320.2,1,329.5,1,360.6,1,382.2,1,435.4,1,464.0,1,468.8,1,403.0,1,351.6,1,252.0,1,188.0,1,146.5,1,152.9,1,148.1,1,165.1,1,177.0,1,206.1,1,244.9,1,228.6,1,253.4,1,241.1,1),dim=c(2,84),dimnames=list(c('Y','D'),1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('Y','D'),1:84))
> 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
Y D M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 96.8 0 1 0 0 0 0 0 0 0 0 0 0 1
2 87.0 0 0 1 0 0 0 0 0 0 0 0 0 2
3 96.3 0 0 0 1 0 0 0 0 0 0 0 0 3
4 107.1 0 0 0 0 1 0 0 0 0 0 0 0 4
5 115.2 0 0 0 0 0 1 0 0 0 0 0 0 5
6 106.1 0 0 0 0 0 0 1 0 0 0 0 0 6
7 89.5 0 0 0 0 0 0 0 1 0 0 0 0 7
8 91.3 0 0 0 0 0 0 0 0 1 0 0 0 8
9 97.6 0 0 0 0 0 0 0 0 0 1 0 0 9
10 100.7 0 0 0 0 0 0 0 0 0 0 1 0 10
11 104.6 0 0 0 0 0 0 0 0 0 0 0 1 11
12 94.7 0 0 0 0 0 0 0 0 0 0 0 0 12
13 101.8 0 1 0 0 0 0 0 0 0 0 0 0 13
14 102.5 0 0 1 0 0 0 0 0 0 0 0 0 14
15 105.3 0 0 0 1 0 0 0 0 0 0 0 0 15
16 110.3 0 0 0 0 1 0 0 0 0 0 0 0 16
17 109.8 0 0 0 0 0 1 0 0 0 0 0 0 17
18 117.3 0 0 0 0 0 0 1 0 0 0 0 0 18
19 118.8 0 0 0 0 0 0 0 1 0 0 0 0 19
20 131.3 0 0 0 0 0 0 0 0 1 0 0 0 20
21 125.9 0 0 0 0 0 0 0 0 0 1 0 0 21
22 133.1 0 0 0 0 0 0 0 0 0 0 1 0 22
23 147.0 0 0 0 0 0 0 0 0 0 0 0 1 23
24 145.8 0 0 0 0 0 0 0 0 0 0 0 0 24
25 164.4 0 1 0 0 0 0 0 0 0 0 0 0 25
26 149.8 0 0 1 0 0 0 0 0 0 0 0 0 26
27 137.7 0 0 0 1 0 0 0 0 0 0 0 0 27
28 151.7 0 0 0 0 1 0 0 0 0 0 0 0 28
29 156.8 0 0 0 0 0 1 0 0 0 0 0 0 29
30 180.0 0 0 0 0 0 0 1 0 0 0 0 0 30
31 180.4 0 0 0 0 0 0 0 1 0 0 0 0 31
32 170.4 0 0 0 0 0 0 0 0 1 0 0 0 32
33 191.6 0 0 0 0 0 0 0 0 0 1 0 0 33
34 199.5 0 0 0 0 0 0 0 0 0 0 1 0 34
35 218.2 0 0 0 0 0 0 0 0 0 0 0 1 35
36 217.5 0 0 0 0 0 0 0 0 0 0 0 0 36
37 205.0 0 1 0 0 0 0 0 0 0 0 0 0 37
38 194.0 0 0 1 0 0 0 0 0 0 0 0 0 38
39 199.3 0 0 0 1 0 0 0 0 0 0 0 0 39
40 219.3 0 0 0 0 1 0 0 0 0 0 0 0 40
41 211.1 0 0 0 0 0 1 0 0 0 0 0 0 41
42 215.2 0 0 0 0 0 0 1 0 0 0 0 0 42
43 240.2 0 0 0 0 0 0 0 1 0 0 0 0 43
44 242.2 0 0 0 0 0 0 0 0 1 0 0 0 44
45 240.7 0 0 0 0 0 0 0 0 0 1 0 0 45
46 255.4 0 0 0 0 0 0 0 0 0 0 1 0 46
47 253.0 0 0 0 0 0 0 0 0 0 0 0 1 47
48 218.2 0 0 0 0 0 0 0 0 0 0 0 0 48
49 203.7 0 1 0 0 0 0 0 0 0 0 0 0 49
50 205.6 0 0 1 0 0 0 0 0 0 0 0 0 50
51 215.6 0 0 0 1 0 0 0 0 0 0 0 0 51
52 188.5 0 0 0 0 1 0 0 0 0 0 0 0 52
53 202.9 0 0 0 0 0 1 0 0 0 0 0 0 53
54 214.0 0 0 0 0 0 0 1 0 0 0 0 0 54
55 230.3 0 0 0 0 0 0 0 1 0 0 0 0 55
56 230.0 0 0 0 0 0 0 0 0 1 0 0 0 56
57 241.0 0 0 0 0 0 0 0 0 0 1 0 0 57
58 259.6 1 0 0 0 0 0 0 0 0 0 1 0 58
59 247.8 1 0 0 0 0 0 0 0 0 0 0 1 59
60 270.3 1 0 0 0 0 0 0 0 0 0 0 0 60
61 289.7 1 1 0 0 0 0 0 0 0 0 0 0 61
62 322.7 1 0 1 0 0 0 0 0 0 0 0 0 62
63 315.0 1 0 0 1 0 0 0 0 0 0 0 0 63
64 320.2 1 0 0 0 1 0 0 0 0 0 0 0 64
65 329.5 1 0 0 0 0 1 0 0 0 0 0 0 65
66 360.6 1 0 0 0 0 0 1 0 0 0 0 0 66
67 382.2 1 0 0 0 0 0 0 1 0 0 0 0 67
68 435.4 1 0 0 0 0 0 0 0 1 0 0 0 68
69 464.0 1 0 0 0 0 0 0 0 0 1 0 0 69
70 468.8 1 0 0 0 0 0 0 0 0 0 1 0 70
71 403.0 1 0 0 0 0 0 0 0 0 0 0 1 71
72 351.6 1 0 0 0 0 0 0 0 0 0 0 0 72
73 252.0 1 1 0 0 0 0 0 0 0 0 0 0 73
74 188.0 1 0 1 0 0 0 0 0 0 0 0 0 74
75 146.5 1 0 0 1 0 0 0 0 0 0 0 0 75
76 152.9 1 0 0 0 1 0 0 0 0 0 0 0 76
77 148.1 1 0 0 0 0 1 0 0 0 0 0 0 77
78 165.1 1 0 0 0 0 0 1 0 0 0 0 0 78
79 177.0 1 0 0 0 0 0 0 1 0 0 0 0 79
80 206.1 1 0 0 0 0 0 0 0 1 0 0 0 80
81 244.9 1 0 0 0 0 0 0 0 0 1 0 0 81
82 228.6 1 0 0 0 0 0 0 0 0 0 1 0 82
83 253.4 1 0 0 0 0 0 0 0 0 0 0 1 83
84 241.1 1 0 0 0 0 0 0 0 0 0 0 0 84
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) D M1 M2 M3 M4
106.454 25.966 -5.103 -16.349 -23.323 -20.554
M5 M6 M7 M8 M9 M10
-19.343 -9.346 -2.891 7.592 19.603 19.477
M11 t
14.674 2.131
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-129.0886 -25.8440 -0.4908 25.7367 167.7109
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 106.4545 29.4790 3.611 0.000569 ***
D 25.9658 26.0368 0.997 0.322067
M1 -5.1033 34.7781 -0.147 0.883760
M2 -16.3489 34.7455 -0.471 0.639437
M3 -23.3231 34.7200 -0.672 0.503956
M4 -20.5544 34.7019 -0.592 0.555549
M5 -19.3428 34.6910 -0.558 0.578913
M6 -9.3456 34.6873 -0.269 0.788397
M7 -2.8912 34.6910 -0.083 0.933818
M8 7.5918 34.7019 0.219 0.827464
M9 19.6033 34.7200 0.565 0.574143
M10 19.4769 34.6333 0.562 0.575657
M11 14.6742 34.6224 0.424 0.672986
t 2.1313 0.5022 4.244 6.62e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 64.77 on 70 degrees of freedom
Multiple R-squared: 0.5483, Adjusted R-squared: 0.4645
F-statistic: 6.537 on 13 and 70 DF, p-value: 6.066e-08
> 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,] 9.478627e-04 1.895725e-03 0.9990521
[2,] 8.468820e-05 1.693764e-04 0.9999153
[3,] 6.967396e-05 1.393479e-04 0.9999303
[4,] 6.992108e-05 1.398422e-04 0.9999301
[5,] 1.685415e-05 3.370831e-05 0.9999831
[6,] 4.792115e-06 9.584230e-06 0.9999952
[7,] 2.446958e-06 4.893915e-06 0.9999976
[8,] 1.919299e-06 3.838599e-06 0.9999981
[9,] 1.504168e-06 3.008337e-06 0.9999985
[10,] 4.281914e-07 8.563829e-07 0.9999996
[11,] 9.163108e-08 1.832622e-07 0.9999999
[12,] 1.803129e-08 3.606258e-08 1.0000000
[13,] 3.492012e-09 6.984023e-09 1.0000000
[14,] 2.168044e-09 4.336088e-09 1.0000000
[15,] 2.227299e-09 4.454597e-09 1.0000000
[16,] 8.525365e-10 1.705073e-09 1.0000000
[17,] 1.297939e-09 2.595878e-09 1.0000000
[18,] 1.507696e-09 3.015392e-09 1.0000000
[19,] 2.469598e-09 4.939197e-09 1.0000000
[20,] 3.836686e-09 7.673372e-09 1.0000000
[21,] 1.064618e-09 2.129236e-09 1.0000000
[22,] 2.837871e-10 5.675743e-10 1.0000000
[23,] 7.301239e-11 1.460248e-10 1.0000000
[24,] 2.264414e-11 4.528829e-11 1.0000000
[25,] 5.126915e-12 1.025383e-11 1.0000000
[26,] 1.193927e-12 2.387853e-12 1.0000000
[27,] 7.911391e-13 1.582278e-12 1.0000000
[28,] 5.618996e-13 1.123799e-12 1.0000000
[29,] 3.454800e-13 6.909601e-13 1.0000000
[30,] 1.645322e-13 3.290644e-13 1.0000000
[31,] 3.885415e-14 7.770829e-14 1.0000000
[32,] 1.652964e-14 3.305929e-14 1.0000000
[33,] 3.774239e-14 7.548477e-14 1.0000000
[34,] 2.286435e-14 4.572870e-14 1.0000000
[35,] 8.664446e-15 1.732889e-14 1.0000000
[36,] 4.724210e-14 9.448421e-14 1.0000000
[37,] 4.805099e-14 9.610198e-14 1.0000000
[38,] 2.884854e-14 5.769709e-14 1.0000000
[39,] 9.414400e-15 1.882880e-14 1.0000000
[40,] 2.668296e-15 5.336592e-15 1.0000000
[41,] 5.934258e-16 1.186852e-15 1.0000000
[42,] 1.912896e-14 3.825793e-14 1.0000000
[43,] 6.082365e-11 1.216473e-10 1.0000000
[44,] 2.032605e-04 4.065210e-04 0.9997967
[45,] 3.272431e-02 6.544863e-02 0.9672757
[46,] 7.324170e-02 1.464834e-01 0.9267583
[47,] 5.942172e-02 1.188434e-01 0.9405783
[48,] 4.543018e-02 9.086035e-02 0.9545698
[49,] 2.872826e-02 5.745652e-02 0.9712717
[50,] 1.978039e-02 3.956078e-02 0.9802196
[51,] 1.439932e-02 2.879864e-02 0.9856007
> postscript(file="/var/www/html/rcomp/tmp/1niq01260703176.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/2alul1260703176.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/3v0wp1260703176.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/4oyuy1260703176.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/5u73o1260703176.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 = 84
Frequency = 1
1 2 3 4 5 6
-6.682498 -7.368213 6.774644 12.674644 17.431787 -3.796784
7 8 9 10 11 12
-28.982498 -39.796784 -47.639641 -46.544527 -39.973098 -37.330241
13 14 15 16 17 18
-27.258256 -17.443970 -9.801113 -9.701113 -13.543970 -18.172542
19 20 21 22 23 24
-25.258256 -25.372542 -44.915399 -39.720284 -23.148856 -11.805999
25 26 27 28 29 30
9.765986 4.280272 -2.976871 6.123129 7.880272 18.951701
31 32 33 34 35 36
10.765986 -11.848299 -4.791156 1.103958 22.475387 34.318244
37 38 39 40 41 42
24.790229 22.904515 33.047372 48.147372 36.604515 28.575943
43 44 45 46 47 48
44.990229 34.375943 18.733086 31.428200 31.699629 9.442486
49 50 51 52 53 54
-2.085529 8.928757 23.771614 -8.228386 2.828757 1.800186
55 56 57 58 59 60
9.514471 -3.399814 -6.542672 -15.913358 -25.041929 10.000928
61 62 63 64 65 66
32.372913 74.487199 71.630056 71.930056 77.887199 96.858627
67 68 69 70 71 72
109.872913 150.458627 164.915770 167.710884 104.582313 65.725170
73 74 75 76 77 78
-30.902845 -85.788559 -122.445702 -120.945702 -129.088559 -124.217130
79 80 81 82 83 84
-120.902845 -104.417130 -79.759988 -98.064873 -70.593445 -70.350588
> postscript(file="/var/www/html/rcomp/tmp/6dtsj1260703176.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 = 84
Frequency = 1
lag(myerror, k = 1) myerror
0 -6.682498 NA
1 -7.368213 -6.682498
2 6.774644 -7.368213
3 12.674644 6.774644
4 17.431787 12.674644
5 -3.796784 17.431787
6 -28.982498 -3.796784
7 -39.796784 -28.982498
8 -47.639641 -39.796784
9 -46.544527 -47.639641
10 -39.973098 -46.544527
11 -37.330241 -39.973098
12 -27.258256 -37.330241
13 -17.443970 -27.258256
14 -9.801113 -17.443970
15 -9.701113 -9.801113
16 -13.543970 -9.701113
17 -18.172542 -13.543970
18 -25.258256 -18.172542
19 -25.372542 -25.258256
20 -44.915399 -25.372542
21 -39.720284 -44.915399
22 -23.148856 -39.720284
23 -11.805999 -23.148856
24 9.765986 -11.805999
25 4.280272 9.765986
26 -2.976871 4.280272
27 6.123129 -2.976871
28 7.880272 6.123129
29 18.951701 7.880272
30 10.765986 18.951701
31 -11.848299 10.765986
32 -4.791156 -11.848299
33 1.103958 -4.791156
34 22.475387 1.103958
35 34.318244 22.475387
36 24.790229 34.318244
37 22.904515 24.790229
38 33.047372 22.904515
39 48.147372 33.047372
40 36.604515 48.147372
41 28.575943 36.604515
42 44.990229 28.575943
43 34.375943 44.990229
44 18.733086 34.375943
45 31.428200 18.733086
46 31.699629 31.428200
47 9.442486 31.699629
48 -2.085529 9.442486
49 8.928757 -2.085529
50 23.771614 8.928757
51 -8.228386 23.771614
52 2.828757 -8.228386
53 1.800186 2.828757
54 9.514471 1.800186
55 -3.399814 9.514471
56 -6.542672 -3.399814
57 -15.913358 -6.542672
58 -25.041929 -15.913358
59 10.000928 -25.041929
60 32.372913 10.000928
61 74.487199 32.372913
62 71.630056 74.487199
63 71.930056 71.630056
64 77.887199 71.930056
65 96.858627 77.887199
66 109.872913 96.858627
67 150.458627 109.872913
68 164.915770 150.458627
69 167.710884 164.915770
70 104.582313 167.710884
71 65.725170 104.582313
72 -30.902845 65.725170
73 -85.788559 -30.902845
74 -122.445702 -85.788559
75 -120.945702 -122.445702
76 -129.088559 -120.945702
77 -124.217130 -129.088559
78 -120.902845 -124.217130
79 -104.417130 -120.902845
80 -79.759988 -104.417130
81 -98.064873 -79.759988
82 -70.593445 -98.064873
83 -70.350588 -70.593445
84 NA -70.350588
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -7.368213 -6.682498
[2,] 6.774644 -7.368213
[3,] 12.674644 6.774644
[4,] 17.431787 12.674644
[5,] -3.796784 17.431787
[6,] -28.982498 -3.796784
[7,] -39.796784 -28.982498
[8,] -47.639641 -39.796784
[9,] -46.544527 -47.639641
[10,] -39.973098 -46.544527
[11,] -37.330241 -39.973098
[12,] -27.258256 -37.330241
[13,] -17.443970 -27.258256
[14,] -9.801113 -17.443970
[15,] -9.701113 -9.801113
[16,] -13.543970 -9.701113
[17,] -18.172542 -13.543970
[18,] -25.258256 -18.172542
[19,] -25.372542 -25.258256
[20,] -44.915399 -25.372542
[21,] -39.720284 -44.915399
[22,] -23.148856 -39.720284
[23,] -11.805999 -23.148856
[24,] 9.765986 -11.805999
[25,] 4.280272 9.765986
[26,] -2.976871 4.280272
[27,] 6.123129 -2.976871
[28,] 7.880272 6.123129
[29,] 18.951701 7.880272
[30,] 10.765986 18.951701
[31,] -11.848299 10.765986
[32,] -4.791156 -11.848299
[33,] 1.103958 -4.791156
[34,] 22.475387 1.103958
[35,] 34.318244 22.475387
[36,] 24.790229 34.318244
[37,] 22.904515 24.790229
[38,] 33.047372 22.904515
[39,] 48.147372 33.047372
[40,] 36.604515 48.147372
[41,] 28.575943 36.604515
[42,] 44.990229 28.575943
[43,] 34.375943 44.990229
[44,] 18.733086 34.375943
[45,] 31.428200 18.733086
[46,] 31.699629 31.428200
[47,] 9.442486 31.699629
[48,] -2.085529 9.442486
[49,] 8.928757 -2.085529
[50,] 23.771614 8.928757
[51,] -8.228386 23.771614
[52,] 2.828757 -8.228386
[53,] 1.800186 2.828757
[54,] 9.514471 1.800186
[55,] -3.399814 9.514471
[56,] -6.542672 -3.399814
[57,] -15.913358 -6.542672
[58,] -25.041929 -15.913358
[59,] 10.000928 -25.041929
[60,] 32.372913 10.000928
[61,] 74.487199 32.372913
[62,] 71.630056 74.487199
[63,] 71.930056 71.630056
[64,] 77.887199 71.930056
[65,] 96.858627 77.887199
[66,] 109.872913 96.858627
[67,] 150.458627 109.872913
[68,] 164.915770 150.458627
[69,] 167.710884 164.915770
[70,] 104.582313 167.710884
[71,] 65.725170 104.582313
[72,] -30.902845 65.725170
[73,] -85.788559 -30.902845
[74,] -122.445702 -85.788559
[75,] -120.945702 -122.445702
[76,] -129.088559 -120.945702
[77,] -124.217130 -129.088559
[78,] -120.902845 -124.217130
[79,] -104.417130 -120.902845
[80,] -79.759988 -104.417130
[81,] -98.064873 -79.759988
[82,] -70.593445 -98.064873
[83,] -70.350588 -70.593445
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -7.368213 -6.682498
2 6.774644 -7.368213
3 12.674644 6.774644
4 17.431787 12.674644
5 -3.796784 17.431787
6 -28.982498 -3.796784
7 -39.796784 -28.982498
8 -47.639641 -39.796784
9 -46.544527 -47.639641
10 -39.973098 -46.544527
11 -37.330241 -39.973098
12 -27.258256 -37.330241
13 -17.443970 -27.258256
14 -9.801113 -17.443970
15 -9.701113 -9.801113
16 -13.543970 -9.701113
17 -18.172542 -13.543970
18 -25.258256 -18.172542
19 -25.372542 -25.258256
20 -44.915399 -25.372542
21 -39.720284 -44.915399
22 -23.148856 -39.720284
23 -11.805999 -23.148856
24 9.765986 -11.805999
25 4.280272 9.765986
26 -2.976871 4.280272
27 6.123129 -2.976871
28 7.880272 6.123129
29 18.951701 7.880272
30 10.765986 18.951701
31 -11.848299 10.765986
32 -4.791156 -11.848299
33 1.103958 -4.791156
34 22.475387 1.103958
35 34.318244 22.475387
36 24.790229 34.318244
37 22.904515 24.790229
38 33.047372 22.904515
39 48.147372 33.047372
40 36.604515 48.147372
41 28.575943 36.604515
42 44.990229 28.575943
43 34.375943 44.990229
44 18.733086 34.375943
45 31.428200 18.733086
46 31.699629 31.428200
47 9.442486 31.699629
48 -2.085529 9.442486
49 8.928757 -2.085529
50 23.771614 8.928757
51 -8.228386 23.771614
52 2.828757 -8.228386
53 1.800186 2.828757
54 9.514471 1.800186
55 -3.399814 9.514471
56 -6.542672 -3.399814
57 -15.913358 -6.542672
58 -25.041929 -15.913358
59 10.000928 -25.041929
60 32.372913 10.000928
61 74.487199 32.372913
62 71.630056 74.487199
63 71.930056 71.630056
64 77.887199 71.930056
65 96.858627 77.887199
66 109.872913 96.858627
67 150.458627 109.872913
68 164.915770 150.458627
69 167.710884 164.915770
70 104.582313 167.710884
71 65.725170 104.582313
72 -30.902845 65.725170
73 -85.788559 -30.902845
74 -122.445702 -85.788559
75 -120.945702 -122.445702
76 -129.088559 -120.945702
77 -124.217130 -129.088559
78 -120.902845 -124.217130
79 -104.417130 -120.902845
80 -79.759988 -104.417130
81 -98.064873 -79.759988
82 -70.593445 -98.064873
83 -70.350588 -70.593445
> 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/74py11260703176.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/81m0b1260703176.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/9a0t21260703176.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/10c2ic1260703176.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/11legy1260703176.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/12f65t1260703176.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/13cuvc1260703176.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/14fnyz1260703176.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/150efv1260703177.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/16ltxm1260703177.tab")
+ }
>
> try(system("convert tmp/1niq01260703176.ps tmp/1niq01260703176.png",intern=TRUE))
character(0)
> try(system("convert tmp/2alul1260703176.ps tmp/2alul1260703176.png",intern=TRUE))
character(0)
> try(system("convert tmp/3v0wp1260703176.ps tmp/3v0wp1260703176.png",intern=TRUE))
character(0)
> try(system("convert tmp/4oyuy1260703176.ps tmp/4oyuy1260703176.png",intern=TRUE))
character(0)
> try(system("convert tmp/5u73o1260703176.ps tmp/5u73o1260703176.png",intern=TRUE))
character(0)
> try(system("convert tmp/6dtsj1260703176.ps tmp/6dtsj1260703176.png",intern=TRUE))
character(0)
> try(system("convert tmp/74py11260703176.ps tmp/74py11260703176.png",intern=TRUE))
character(0)
> try(system("convert tmp/81m0b1260703176.ps tmp/81m0b1260703176.png",intern=TRUE))
character(0)
> try(system("convert tmp/9a0t21260703176.ps tmp/9a0t21260703176.png",intern=TRUE))
character(0)
> try(system("convert tmp/10c2ic1260703176.ps tmp/10c2ic1260703176.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.689 1.571 3.284