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(98.5,0,97.0,0,103.3,0,99.6,0,100.1,0,102.9,0,95.9,0,94.5,0,107.4,0,116.0,0,102.8,0,99.8,0,109.6,0,103.0,0,111.6,0,106.3,0,97.9,0,108.8,0,103.9,0,101.2,0,122.9,0,123.9,0,111.7,0,120.9,0,99.6,0,103.3,0,119.4,0,106.5,0,101.9,0,124.6,0,106.5,0,107.8,0,127.4,0,120.1,0,118.5,0,127.7,0,107.7,0,104.5,0,118.8,0,110.3,0,109.6,0,119.1,0,96.5,0,106.7,0,126.3,0,116.2,0,118.8,0,115.2,0,110.0,0,111.4,0,129.6,0,108.1,0,117.8,0,122.9,0,100.6,0,111.8,0,127.0,0,128.6,0,124.8,0,118.5,0,114.7,0,112.6,0,128.7,0,111.0,0,115.8,0,126.0,0,111.1,1,113.2,1,120.1,1,130.6,1,124.0,1,119.4,1,116.7,1,116.5,1,119.6,1,126.5,1,111.3,1,123.5,1,114.2,1,103.7,1,129.5,1),dim=c(2,81),dimnames=list(c('Y','X'),1:81))
> y <- array(NA,dim=c(2,81),dimnames=list(c('Y','X'),1:81))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 98.5 0 1 0 0 0 0 0 0 0 0 0 0
2 97.0 0 0 1 0 0 0 0 0 0 0 0 0
3 103.3 0 0 0 1 0 0 0 0 0 0 0 0
4 99.6 0 0 0 0 1 0 0 0 0 0 0 0
5 100.1 0 0 0 0 0 1 0 0 0 0 0 0
6 102.9 0 0 0 0 0 0 1 0 0 0 0 0
7 95.9 0 0 0 0 0 0 0 1 0 0 0 0
8 94.5 0 0 0 0 0 0 0 0 1 0 0 0
9 107.4 0 0 0 0 0 0 0 0 0 1 0 0
10 116.0 0 0 0 0 0 0 0 0 0 0 1 0
11 102.8 0 0 0 0 0 0 0 0 0 0 0 1
12 99.8 0 0 0 0 0 0 0 0 0 0 0 0
13 109.6 0 1 0 0 0 0 0 0 0 0 0 0
14 103.0 0 0 1 0 0 0 0 0 0 0 0 0
15 111.6 0 0 0 1 0 0 0 0 0 0 0 0
16 106.3 0 0 0 0 1 0 0 0 0 0 0 0
17 97.9 0 0 0 0 0 1 0 0 0 0 0 0
18 108.8 0 0 0 0 0 0 1 0 0 0 0 0
19 103.9 0 0 0 0 0 0 0 1 0 0 0 0
20 101.2 0 0 0 0 0 0 0 0 1 0 0 0
21 122.9 0 0 0 0 0 0 0 0 0 1 0 0
22 123.9 0 0 0 0 0 0 0 0 0 0 1 0
23 111.7 0 0 0 0 0 0 0 0 0 0 0 1
24 120.9 0 0 0 0 0 0 0 0 0 0 0 0
25 99.6 0 1 0 0 0 0 0 0 0 0 0 0
26 103.3 0 0 1 0 0 0 0 0 0 0 0 0
27 119.4 0 0 0 1 0 0 0 0 0 0 0 0
28 106.5 0 0 0 0 1 0 0 0 0 0 0 0
29 101.9 0 0 0 0 0 1 0 0 0 0 0 0
30 124.6 0 0 0 0 0 0 1 0 0 0 0 0
31 106.5 0 0 0 0 0 0 0 1 0 0 0 0
32 107.8 0 0 0 0 0 0 0 0 1 0 0 0
33 127.4 0 0 0 0 0 0 0 0 0 1 0 0
34 120.1 0 0 0 0 0 0 0 0 0 0 1 0
35 118.5 0 0 0 0 0 0 0 0 0 0 0 1
36 127.7 0 0 0 0 0 0 0 0 0 0 0 0
37 107.7 0 1 0 0 0 0 0 0 0 0 0 0
38 104.5 0 0 1 0 0 0 0 0 0 0 0 0
39 118.8 0 0 0 1 0 0 0 0 0 0 0 0
40 110.3 0 0 0 0 1 0 0 0 0 0 0 0
41 109.6 0 0 0 0 0 1 0 0 0 0 0 0
42 119.1 0 0 0 0 0 0 1 0 0 0 0 0
43 96.5 0 0 0 0 0 0 0 1 0 0 0 0
44 106.7 0 0 0 0 0 0 0 0 1 0 0 0
45 126.3 0 0 0 0 0 0 0 0 0 1 0 0
46 116.2 0 0 0 0 0 0 0 0 0 0 1 0
47 118.8 0 0 0 0 0 0 0 0 0 0 0 1
48 115.2 0 0 0 0 0 0 0 0 0 0 0 0
49 110.0 0 1 0 0 0 0 0 0 0 0 0 0
50 111.4 0 0 1 0 0 0 0 0 0 0 0 0
51 129.6 0 0 0 1 0 0 0 0 0 0 0 0
52 108.1 0 0 0 0 1 0 0 0 0 0 0 0
53 117.8 0 0 0 0 0 1 0 0 0 0 0 0
54 122.9 0 0 0 0 0 0 1 0 0 0 0 0
55 100.6 0 0 0 0 0 0 0 1 0 0 0 0
56 111.8 0 0 0 0 0 0 0 0 1 0 0 0
57 127.0 0 0 0 0 0 0 0 0 0 1 0 0
58 128.6 0 0 0 0 0 0 0 0 0 0 1 0
59 124.8 0 0 0 0 0 0 0 0 0 0 0 1
60 118.5 0 0 0 0 0 0 0 0 0 0 0 0
61 114.7 0 1 0 0 0 0 0 0 0 0 0 0
62 112.6 0 0 1 0 0 0 0 0 0 0 0 0
63 128.7 0 0 0 1 0 0 0 0 0 0 0 0
64 111.0 0 0 0 0 1 0 0 0 0 0 0 0
65 115.8 0 0 0 0 0 1 0 0 0 0 0 0
66 126.0 0 0 0 0 0 0 1 0 0 0 0 0
67 111.1 1 0 0 0 0 0 0 1 0 0 0 0
68 113.2 1 0 0 0 0 0 0 0 1 0 0 0
69 120.1 1 0 0 0 0 0 0 0 0 1 0 0
70 130.6 1 0 0 0 0 0 0 0 0 0 1 0
71 124.0 1 0 0 0 0 0 0 0 0 0 0 1
72 119.4 1 0 0 0 0 0 0 0 0 0 0 0
73 116.7 1 1 0 0 0 0 0 0 0 0 0 0
74 116.5 1 0 1 0 0 0 0 0 0 0 0 0
75 119.6 1 0 0 1 0 0 0 0 0 0 0 0
76 126.5 1 0 0 0 1 0 0 0 0 0 0 0
77 111.3 1 0 0 0 0 1 0 0 0 0 0 0
78 123.5 1 0 0 0 0 0 1 0 0 0 0 0
79 114.2 1 0 0 0 0 0 0 1 0 0 0 0
80 103.7 1 0 0 0 0 0 0 0 1 0 0 0
81 129.5 1 0 0 0 0 0 0 0 0 1 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
115.674 7.456 -8.625 -9.839 1.975 -6.982
M5 M6 M7 M8 M9 M10
-8.968 1.518 -13.704 -12.247 5.139 5.650
M11
-0.150
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-15.874 -5.124 1.608 4.530 12.026
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 115.674 2.959 39.090 < 2e-16 ***
X 7.456 2.084 3.577 0.000645 ***
M1 -8.625 4.005 -2.153 0.034831 *
M2 -9.839 4.005 -2.457 0.016582 *
M3 1.975 4.005 0.493 0.623493
M4 -6.982 4.005 -1.743 0.085805 .
M5 -8.968 4.005 -2.239 0.028427 *
M6 1.518 4.005 0.379 0.705858
M7 -13.704 4.013 -3.415 0.001078 **
M8 -12.247 4.013 -3.052 0.003239 **
M9 5.139 4.013 1.281 0.204665
M10 5.650 4.156 1.359 0.178483
M11 -0.150 4.156 -0.036 0.971314
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.198 on 68 degrees of freedom
Multiple R-squared: 0.5401, Adjusted R-squared: 0.459
F-statistic: 6.655 on 12 and 68 DF, p-value: 1.017e-07
> 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.7592013 0.48159731 0.24079865
[2,] 0.6661291 0.66774176 0.33387088
[3,] 0.6536458 0.69270848 0.34635424
[4,] 0.6421838 0.71563247 0.35781623
[5,] 0.6128038 0.77439241 0.38719620
[6,] 0.8033666 0.39326690 0.19663345
[7,] 0.7809543 0.43809140 0.21904570
[8,] 0.8004043 0.39919136 0.19959568
[9,] 0.9451321 0.10973571 0.05486786
[10,] 0.9490285 0.10194309 0.05097154
[11,] 0.9363487 0.12730251 0.06365126
[12,] 0.9527207 0.09455852 0.04727926
[13,] 0.9414879 0.11702423 0.05851212
[14,] 0.9462851 0.10742972 0.05371486
[15,] 0.9825236 0.03495277 0.01747639
[16,] 0.9770915 0.04581701 0.02290850
[17,] 0.9750646 0.04987089 0.02493545
[18,] 0.9781956 0.04360878 0.02180439
[19,] 0.9679484 0.06410325 0.03205162
[20,] 0.9679227 0.06415454 0.03207727
[21,] 0.9895119 0.02097617 0.01048809
[22,] 0.9861786 0.02764285 0.01382142
[23,] 0.9849227 0.03015461 0.01507731
[24,] 0.9834894 0.03302114 0.01651057
[25,] 0.9782973 0.04340548 0.02170274
[26,] 0.9761820 0.04763593 0.02381796
[27,] 0.9700134 0.05997318 0.02998659
[28,] 0.9778007 0.04439852 0.02219926
[29,] 0.9677044 0.06459129 0.03229565
[30,] 0.9567954 0.08640918 0.04320459
[31,] 0.9738289 0.05234218 0.02617109
[32,] 0.9684239 0.06315219 0.03157610
[33,] 0.9530979 0.09380425 0.04690212
[34,] 0.9425307 0.11493869 0.05746934
[35,] 0.9286692 0.14266168 0.07133084
[36,] 0.9444601 0.11107982 0.05553991
[37,] 0.9555275 0.08894500 0.04447250
[38,] 0.9577267 0.08454664 0.04227332
[39,] 0.9391912 0.12161763 0.06080881
[40,] 0.9638148 0.07237039 0.03618519
[41,] 0.9501513 0.09969745 0.04984872
[42,] 0.9247675 0.15046502 0.07523251
[43,] 0.8891272 0.22174570 0.11087285
[44,] 0.8463443 0.30731135 0.15365567
[45,] 0.7680480 0.46390407 0.23195204
[46,] 0.6834264 0.63314717 0.31657359
[47,] 0.5902996 0.81940079 0.40970040
[48,] 0.6025155 0.79496899 0.39748450
[49,] 0.8009108 0.39817841 0.19908920
[50,] 0.6604507 0.67909862 0.33954931
> postscript(file="/var/www/html/rcomp/tmp/10zv81229431425.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/2gg7m1229431425.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/3ivdh1229431425.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/409nt1229431425.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/59raw1229431425.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 = 81
Frequency = 1
1 2 3 4 5 6
-8.5491873 -8.8349016 -14.3491873 -9.0920445 -6.6063302 -14.2920445
7 8 9 10 11 12
-6.0698033 -8.9269461 -13.4126604 -5.3240519 -12.7240519 -15.8740519
13 14 15 16 17 18
2.5508127 -2.8349016 -6.0491873 -2.3920445 -8.8063302 -8.3920445
19 20 21 22 23 24
1.9301967 -2.2269461 2.0873396 2.5759481 -3.8240519 5.2259481
25 26 27 28 29 30
-7.4491873 -2.5349016 1.7508127 -2.1920445 -4.8063302 7.4079555
31 32 33 34 35 36
4.5301967 4.3730539 6.5873396 -1.2240519 2.9759481 12.0259481
37 38 39 40 41 42
0.6508127 -1.3349016 1.1508127 1.6079555 2.8936698 1.9079555
43 44 45 46 47 48
-5.4698033 3.2730539 5.4873396 -5.1240519 3.2759481 -0.4740519
49 50 51 52 53 54
2.9508127 5.5650984 11.9508127 -0.5920445 11.0936698 5.7079555
55 56 57 58 59 60
-1.3698033 8.3730539 6.1873396 7.2759481 9.2759481 2.8259481
61 62 63 64 65 66
7.6508127 6.7650984 11.0508127 2.3079555 9.0936698 8.8079555
67 68 69 70 71 72
1.6745081 2.3173653 -8.1683490 1.8202595 1.0202595 -3.7297405
73 74 75 76 77 78
2.1951240 3.2094098 -5.5048760 10.3522669 -2.8620188 -1.1477331
79 80 81
4.7745081 -7.1826347 1.2316510
> postscript(file="/var/www/html/rcomp/tmp/6x6io1229431425.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 = 81
Frequency = 1
lag(myerror, k = 1) myerror
0 -8.5491873 NA
1 -8.8349016 -8.5491873
2 -14.3491873 -8.8349016
3 -9.0920445 -14.3491873
4 -6.6063302 -9.0920445
5 -14.2920445 -6.6063302
6 -6.0698033 -14.2920445
7 -8.9269461 -6.0698033
8 -13.4126604 -8.9269461
9 -5.3240519 -13.4126604
10 -12.7240519 -5.3240519
11 -15.8740519 -12.7240519
12 2.5508127 -15.8740519
13 -2.8349016 2.5508127
14 -6.0491873 -2.8349016
15 -2.3920445 -6.0491873
16 -8.8063302 -2.3920445
17 -8.3920445 -8.8063302
18 1.9301967 -8.3920445
19 -2.2269461 1.9301967
20 2.0873396 -2.2269461
21 2.5759481 2.0873396
22 -3.8240519 2.5759481
23 5.2259481 -3.8240519
24 -7.4491873 5.2259481
25 -2.5349016 -7.4491873
26 1.7508127 -2.5349016
27 -2.1920445 1.7508127
28 -4.8063302 -2.1920445
29 7.4079555 -4.8063302
30 4.5301967 7.4079555
31 4.3730539 4.5301967
32 6.5873396 4.3730539
33 -1.2240519 6.5873396
34 2.9759481 -1.2240519
35 12.0259481 2.9759481
36 0.6508127 12.0259481
37 -1.3349016 0.6508127
38 1.1508127 -1.3349016
39 1.6079555 1.1508127
40 2.8936698 1.6079555
41 1.9079555 2.8936698
42 -5.4698033 1.9079555
43 3.2730539 -5.4698033
44 5.4873396 3.2730539
45 -5.1240519 5.4873396
46 3.2759481 -5.1240519
47 -0.4740519 3.2759481
48 2.9508127 -0.4740519
49 5.5650984 2.9508127
50 11.9508127 5.5650984
51 -0.5920445 11.9508127
52 11.0936698 -0.5920445
53 5.7079555 11.0936698
54 -1.3698033 5.7079555
55 8.3730539 -1.3698033
56 6.1873396 8.3730539
57 7.2759481 6.1873396
58 9.2759481 7.2759481
59 2.8259481 9.2759481
60 7.6508127 2.8259481
61 6.7650984 7.6508127
62 11.0508127 6.7650984
63 2.3079555 11.0508127
64 9.0936698 2.3079555
65 8.8079555 9.0936698
66 1.6745081 8.8079555
67 2.3173653 1.6745081
68 -8.1683490 2.3173653
69 1.8202595 -8.1683490
70 1.0202595 1.8202595
71 -3.7297405 1.0202595
72 2.1951240 -3.7297405
73 3.2094098 2.1951240
74 -5.5048760 3.2094098
75 10.3522669 -5.5048760
76 -2.8620188 10.3522669
77 -1.1477331 -2.8620188
78 4.7745081 -1.1477331
79 -7.1826347 4.7745081
80 1.2316510 -7.1826347
81 NA 1.2316510
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -8.8349016 -8.5491873
[2,] -14.3491873 -8.8349016
[3,] -9.0920445 -14.3491873
[4,] -6.6063302 -9.0920445
[5,] -14.2920445 -6.6063302
[6,] -6.0698033 -14.2920445
[7,] -8.9269461 -6.0698033
[8,] -13.4126604 -8.9269461
[9,] -5.3240519 -13.4126604
[10,] -12.7240519 -5.3240519
[11,] -15.8740519 -12.7240519
[12,] 2.5508127 -15.8740519
[13,] -2.8349016 2.5508127
[14,] -6.0491873 -2.8349016
[15,] -2.3920445 -6.0491873
[16,] -8.8063302 -2.3920445
[17,] -8.3920445 -8.8063302
[18,] 1.9301967 -8.3920445
[19,] -2.2269461 1.9301967
[20,] 2.0873396 -2.2269461
[21,] 2.5759481 2.0873396
[22,] -3.8240519 2.5759481
[23,] 5.2259481 -3.8240519
[24,] -7.4491873 5.2259481
[25,] -2.5349016 -7.4491873
[26,] 1.7508127 -2.5349016
[27,] -2.1920445 1.7508127
[28,] -4.8063302 -2.1920445
[29,] 7.4079555 -4.8063302
[30,] 4.5301967 7.4079555
[31,] 4.3730539 4.5301967
[32,] 6.5873396 4.3730539
[33,] -1.2240519 6.5873396
[34,] 2.9759481 -1.2240519
[35,] 12.0259481 2.9759481
[36,] 0.6508127 12.0259481
[37,] -1.3349016 0.6508127
[38,] 1.1508127 -1.3349016
[39,] 1.6079555 1.1508127
[40,] 2.8936698 1.6079555
[41,] 1.9079555 2.8936698
[42,] -5.4698033 1.9079555
[43,] 3.2730539 -5.4698033
[44,] 5.4873396 3.2730539
[45,] -5.1240519 5.4873396
[46,] 3.2759481 -5.1240519
[47,] -0.4740519 3.2759481
[48,] 2.9508127 -0.4740519
[49,] 5.5650984 2.9508127
[50,] 11.9508127 5.5650984
[51,] -0.5920445 11.9508127
[52,] 11.0936698 -0.5920445
[53,] 5.7079555 11.0936698
[54,] -1.3698033 5.7079555
[55,] 8.3730539 -1.3698033
[56,] 6.1873396 8.3730539
[57,] 7.2759481 6.1873396
[58,] 9.2759481 7.2759481
[59,] 2.8259481 9.2759481
[60,] 7.6508127 2.8259481
[61,] 6.7650984 7.6508127
[62,] 11.0508127 6.7650984
[63,] 2.3079555 11.0508127
[64,] 9.0936698 2.3079555
[65,] 8.8079555 9.0936698
[66,] 1.6745081 8.8079555
[67,] 2.3173653 1.6745081
[68,] -8.1683490 2.3173653
[69,] 1.8202595 -8.1683490
[70,] 1.0202595 1.8202595
[71,] -3.7297405 1.0202595
[72,] 2.1951240 -3.7297405
[73,] 3.2094098 2.1951240
[74,] -5.5048760 3.2094098
[75,] 10.3522669 -5.5048760
[76,] -2.8620188 10.3522669
[77,] -1.1477331 -2.8620188
[78,] 4.7745081 -1.1477331
[79,] -7.1826347 4.7745081
[80,] 1.2316510 -7.1826347
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -8.8349016 -8.5491873
2 -14.3491873 -8.8349016
3 -9.0920445 -14.3491873
4 -6.6063302 -9.0920445
5 -14.2920445 -6.6063302
6 -6.0698033 -14.2920445
7 -8.9269461 -6.0698033
8 -13.4126604 -8.9269461
9 -5.3240519 -13.4126604
10 -12.7240519 -5.3240519
11 -15.8740519 -12.7240519
12 2.5508127 -15.8740519
13 -2.8349016 2.5508127
14 -6.0491873 -2.8349016
15 -2.3920445 -6.0491873
16 -8.8063302 -2.3920445
17 -8.3920445 -8.8063302
18 1.9301967 -8.3920445
19 -2.2269461 1.9301967
20 2.0873396 -2.2269461
21 2.5759481 2.0873396
22 -3.8240519 2.5759481
23 5.2259481 -3.8240519
24 -7.4491873 5.2259481
25 -2.5349016 -7.4491873
26 1.7508127 -2.5349016
27 -2.1920445 1.7508127
28 -4.8063302 -2.1920445
29 7.4079555 -4.8063302
30 4.5301967 7.4079555
31 4.3730539 4.5301967
32 6.5873396 4.3730539
33 -1.2240519 6.5873396
34 2.9759481 -1.2240519
35 12.0259481 2.9759481
36 0.6508127 12.0259481
37 -1.3349016 0.6508127
38 1.1508127 -1.3349016
39 1.6079555 1.1508127
40 2.8936698 1.6079555
41 1.9079555 2.8936698
42 -5.4698033 1.9079555
43 3.2730539 -5.4698033
44 5.4873396 3.2730539
45 -5.1240519 5.4873396
46 3.2759481 -5.1240519
47 -0.4740519 3.2759481
48 2.9508127 -0.4740519
49 5.5650984 2.9508127
50 11.9508127 5.5650984
51 -0.5920445 11.9508127
52 11.0936698 -0.5920445
53 5.7079555 11.0936698
54 -1.3698033 5.7079555
55 8.3730539 -1.3698033
56 6.1873396 8.3730539
57 7.2759481 6.1873396
58 9.2759481 7.2759481
59 2.8259481 9.2759481
60 7.6508127 2.8259481
61 6.7650984 7.6508127
62 11.0508127 6.7650984
63 2.3079555 11.0508127
64 9.0936698 2.3079555
65 8.8079555 9.0936698
66 1.6745081 8.8079555
67 2.3173653 1.6745081
68 -8.1683490 2.3173653
69 1.8202595 -8.1683490
70 1.0202595 1.8202595
71 -3.7297405 1.0202595
72 2.1951240 -3.7297405
73 3.2094098 2.1951240
74 -5.5048760 3.2094098
75 10.3522669 -5.5048760
76 -2.8620188 10.3522669
77 -1.1477331 -2.8620188
78 4.7745081 -1.1477331
79 -7.1826347 4.7745081
80 1.2316510 -7.1826347
> 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/7altr1229431425.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/8dl1a1229431425.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/9t9ha1229431425.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/10kjrc1229431425.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/11peov1229431425.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/120v1x1229431425.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/13un8b1229431425.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/1455dc1229431425.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/15degm1229431425.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/16jj341229431426.tab")
+ }
>
> system("convert tmp/10zv81229431425.ps tmp/10zv81229431425.png")
> system("convert tmp/2gg7m1229431425.ps tmp/2gg7m1229431425.png")
> system("convert tmp/3ivdh1229431425.ps tmp/3ivdh1229431425.png")
> system("convert tmp/409nt1229431425.ps tmp/409nt1229431425.png")
> system("convert tmp/59raw1229431425.ps tmp/59raw1229431425.png")
> system("convert tmp/6x6io1229431425.ps tmp/6x6io1229431425.png")
> system("convert tmp/7altr1229431425.ps tmp/7altr1229431425.png")
> system("convert tmp/8dl1a1229431425.ps tmp/8dl1a1229431425.png")
> system("convert tmp/9t9ha1229431425.ps tmp/9t9ha1229431425.png")
> system("convert tmp/10kjrc1229431425.ps tmp/10kjrc1229431425.png")
>
>
> proc.time()
user system elapsed
2.723 1.573 3.269