R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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(1,1,4,1,1,0,0,1,4,0,0,0,1,1,0,1,1,0,1,1,0,0,1,0,0,1,4,1,1,1,0,0,4,0,1,0,0,1,2,0,1,0,0,0,0,1,1,0,1,1,1,1,1,0,0,1,0,0,0,2,1,1,2,1,1,1,0,0,2,1,0,0,1,1,3,1,0,0,1,1,0,0,0,0,0,0,1,1,1,0,1,0,0,1,1,4,0,0,0,0,0,1,0,0,0,1,1,0,1,1,4,0,1,1,0,1,0,1,1,4,1,1,0,1,1,4,1,1,0,1,1,0,0,0,0,0,1,4,0,1,0,1,1,0,1,1,4,0,0,4,0,1,0,1,1,1,0,1,0,0,0,4,0,1,0,0,1,2,0,1,0,0,1,4,0,0,4,0,0,0,0,1,0,1,1,4,1,1,0,1,0,0,0,0,2,0,1,0,0,1,0,0,0,0,1,1,4,1,1,4,0,1,2,0,1,0,0,1,0,0,1,4,1,1,0,1,0,0,0,0,1,1,1,2,1,0,0,1,1,2,0,0,0,0,0,4,0,0,4,1,0,0,0,0,0,0,0,4,1,0,0,1,1,4,0,0,2,0,0,2,1,1,0,1,1,0,1,1,4,0,1,0,1,1,0,1,1,0,1,1,4,1,1,4,0,0,0,0,0,0,1,1,2,0,0,1,0,0,0,0,0,2,0,1,1),dim=c(3,105),dimnames=list(c('pre','post1','post2'),1:105))
> y <- array(NA,dim=c(3,105),dimnames=list(c('pre','post1','post2'),1:105))
> 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 = 'Do not include Seasonal Dummies'
> par1 = '3'
> library(lattice)
> library(lmtest)
Loading required package: zoo
> 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
post2 pre post1
1 4 1 1
2 0 1 1
3 4 0 1
4 0 0 0
5 0 1 1
6 0 1 1
7 0 1 1
8 0 0 1
9 4 0 1
10 1 1 1
11 4 0 0
12 0 0 1
13 2 0 1
14 0 0 1
15 0 0 0
16 0 1 1
17 1 1 1
18 0 1 1
19 0 0 1
20 2 0 0
21 2 1 1
22 1 1 1
23 2 0 0
24 0 1 0
25 3 1 1
26 0 1 0
27 0 1 1
28 0 0 0
29 1 0 0
30 0 1 1
31 0 1 0
32 4 1 1
33 0 0 0
34 1 0 0
35 0 0 0
36 0 1 1
37 4 1 1
38 1 0 1
39 0 0 1
40 4 1 1
41 0 1 1
42 4 1 1
43 0 1 1
44 0 1 1
45 0 0 0
46 4 0 1
47 0 0 1
48 0 1 1
49 4 1 1
50 4 0 0
51 0 0 1
52 1 1 1
53 0 0 1
54 4 0 0
55 0 0 1
56 2 0 1
57 0 0 1
58 4 0 1
59 4 0 0
60 0 0 0
61 0 0 1
62 4 1 1
63 0 1 1
64 0 1 0
65 2 0 0
66 0 0 1
67 0 0 1
68 0 0 0
69 4 1 1
70 4 1 1
71 2 0 1
72 0 0 1
73 0 0 1
74 4 0 1
75 0 1 1
76 0 1 0
77 1 0 0
78 2 1 1
79 0 1 0
80 2 1 1
81 0 0 0
82 4 0 0
83 4 0 0
84 0 1 0
85 0 0 0
86 4 0 0
87 0 1 0
88 4 1 1
89 2 0 0
90 2 0 0
91 0 1 1
92 0 1 1
93 4 1 1
94 0 0 1
95 0 1 1
96 0 1 1
97 4 1 1
98 4 1 1
99 0 0 0
100 0 0 0
101 2 1 1
102 1 0 0
103 0 0 0
104 2 0 0
105 1 0 1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) pre post1
1.16166 -0.01788 0.26719
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.4288 -1.4110 -1.1438 0.8383 2.8383
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.16166 0.28142 4.128 7.5e-05 ***
pre -0.01788 0.35477 -0.050 0.96
post1 0.26719 0.36831 0.725 0.47
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.673 on 102 degrees of freedom
Multiple R-squared: 0.005737, Adjusted R-squared: -0.01376
F-statistic: 0.2943 on 2 and 102 DF, p-value: 0.7457
> 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.7674029 0.4651943 0.23259713
[2,] 0.6659365 0.6681270 0.33406351
[3,] 0.8031891 0.3936217 0.19681087
[4,] 0.7902668 0.4194664 0.20973321
[5,] 0.6997848 0.6004304 0.30021521
[6,] 0.8180850 0.3638301 0.18191503
[7,] 0.8451504 0.3096993 0.15484964
[8,] 0.7831916 0.4336168 0.21680840
[9,] 0.7802894 0.4394212 0.21971062
[10,] 0.7641290 0.4717420 0.23587102
[11,] 0.7136279 0.5727442 0.28637209
[12,] 0.6418343 0.7163315 0.35816573
[13,] 0.5851471 0.8297059 0.41485295
[14,] 0.5692810 0.8614380 0.43071902
[15,] 0.5033027 0.9933946 0.49669732
[16,] 0.4636158 0.9272317 0.53638417
[17,] 0.3926991 0.7853982 0.60730090
[18,] 0.3316417 0.6632834 0.66835828
[19,] 0.2875596 0.5751192 0.71244042
[20,] 0.3207197 0.6414394 0.67928030
[21,] 0.2767802 0.5535605 0.72321975
[22,] 0.2443463 0.4886925 0.75565374
[23,] 0.2224342 0.4448684 0.77756581
[24,] 0.1759271 0.3518541 0.82407293
[25,] 0.1526277 0.3052555 0.84737227
[26,] 0.1246124 0.2492249 0.87538756
[27,] 0.2203634 0.4407268 0.77963661
[28,] 0.1946621 0.3893241 0.80533793
[29,] 0.1541416 0.3082832 0.84585839
[30,] 0.1330093 0.2660187 0.86699066
[31,] 0.1183496 0.2366991 0.88165044
[32,] 0.1905901 0.3811801 0.80940994
[33,] 0.1553954 0.3107908 0.84460459
[34,] 0.1476400 0.2952799 0.85236003
[35,] 0.2146948 0.4293895 0.78530524
[36,] 0.1995590 0.3991179 0.80044103
[37,] 0.2691011 0.5382021 0.73089894
[38,] 0.2534540 0.5069080 0.74654601
[39,] 0.2384063 0.4768126 0.76159369
[40,] 0.2112592 0.4225185 0.78874076
[41,] 0.2699379 0.5398758 0.73006208
[42,] 0.2597301 0.5194602 0.74026990
[43,] 0.2457857 0.4915714 0.75421429
[44,] 0.3146678 0.6293356 0.68533218
[45,] 0.4231548 0.8463096 0.57684518
[46,] 0.4082750 0.8165501 0.59172495
[47,] 0.3578652 0.7157304 0.64213481
[48,] 0.3431820 0.6863639 0.65681804
[49,] 0.4449770 0.8899540 0.55502299
[50,] 0.4299102 0.8598204 0.57008980
[51,] 0.3802870 0.7605741 0.61971295
[52,] 0.3675642 0.7351284 0.63243579
[53,] 0.4340072 0.8680145 0.56599277
[54,] 0.5387404 0.9225193 0.46125965
[55,] 0.5071611 0.9856779 0.49283894
[56,] 0.4934819 0.9869638 0.50651812
[57,] 0.5612353 0.8775293 0.43876466
[58,] 0.5501007 0.8997986 0.44989930
[59,] 0.5202714 0.9594573 0.47972864
[60,] 0.4766896 0.9533793 0.52331035
[61,] 0.4648233 0.9296467 0.53517666
[62,] 0.4598343 0.9196686 0.54016572
[63,] 0.4285364 0.8570727 0.57146364
[64,] 0.4910321 0.9820641 0.50896793
[65,] 0.5605305 0.8789391 0.43946953
[66,] 0.5029425 0.9941150 0.49705749
[67,] 0.5062092 0.9875816 0.49379081
[68,] 0.5308184 0.9383633 0.46918164
[69,] 0.5615805 0.8768390 0.43841948
[70,] 0.5524928 0.8950143 0.44750716
[71,] 0.5089766 0.9820467 0.49102337
[72,] 0.4452757 0.8905515 0.55472426
[73,] 0.3845328 0.7690656 0.61546721
[74,] 0.3480367 0.6960733 0.65196333
[75,] 0.2908585 0.5817170 0.70914148
[76,] 0.2669537 0.5339073 0.73304635
[77,] 0.3534907 0.7069814 0.64650930
[78,] 0.4774581 0.9549162 0.52254190
[79,] 0.4456076 0.8912151 0.55439243
[80,] 0.4031058 0.8062116 0.59689421
[81,] 0.5603907 0.8792185 0.43960925
[82,] 0.5721084 0.8557831 0.42789157
[83,] 0.6331378 0.7337244 0.36686222
[84,] 0.5832539 0.8334922 0.41674612
[85,] 0.5429235 0.9141530 0.45707650
[86,] 0.5599180 0.8801640 0.44008202
[87,] 0.6182637 0.7634727 0.38173634
[88,] 0.6551321 0.6897357 0.34486785
[89,] 0.5756349 0.8487302 0.42436509
[90,] 0.6657098 0.6685805 0.33429023
[91,] 0.9125432 0.1749137 0.08745683
[92,] 0.8756311 0.2487377 0.12436885
[93,] 0.8955571 0.2088858 0.10444290
[94,] 0.8160357 0.3679287 0.18396435
> postscript(file="/var/wessaorg/rcomp/tmp/1ro1k1323987841.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/2rewj1323987841.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/383vs1323987841.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/4gni81323987841.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/5chkv1323987841.ps",horizontal=F,onefile=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 = 105
Frequency = 1
1 2 3 4 5 6 7
2.5890295 -1.4109705 2.5711458 -1.1616597 -1.4109705 -1.4109705 -1.4109705
8 9 10 11 12 13 14
-1.4288542 2.5711458 -0.4109705 2.8383403 -1.4288542 0.5711458 -1.4288542
15 16 17 18 19 20 21
-1.1616597 -1.4109705 -0.4109705 -1.4109705 -1.4288542 0.8383403 0.5890295
22 23 24 25 26 27 28
-0.4109705 0.8383403 -1.1437761 1.5890295 -1.1437761 -1.4109705 -1.1616597
29 30 31 32 33 34 35
-0.1616597 -1.4109705 -1.1437761 2.5890295 -1.1616597 -0.1616597 -1.1616597
36 37 38 39 40 41 42
-1.4109705 2.5890295 -0.4288542 -1.4288542 2.5890295 -1.4109705 2.5890295
43 44 45 46 47 48 49
-1.4109705 -1.4109705 -1.1616597 2.5711458 -1.4288542 -1.4109705 2.5890295
50 51 52 53 54 55 56
2.8383403 -1.4288542 -0.4109705 -1.4288542 2.8383403 -1.4288542 0.5711458
57 58 59 60 61 62 63
-1.4288542 2.5711458 2.8383403 -1.1616597 -1.4288542 2.5890295 -1.4109705
64 65 66 67 68 69 70
-1.1437761 0.8383403 -1.4288542 -1.4288542 -1.1616597 2.5890295 2.5890295
71 72 73 74 75 76 77
0.5711458 -1.4288542 -1.4288542 2.5711458 -1.4109705 -1.1437761 -0.1616597
78 79 80 81 82 83 84
0.5890295 -1.1437761 0.5890295 -1.1616597 2.8383403 2.8383403 -1.1437761
85 86 87 88 89 90 91
-1.1616597 2.8383403 -1.1437761 2.5890295 0.8383403 0.8383403 -1.4109705
92 93 94 95 96 97 98
-1.4109705 2.5890295 -1.4288542 -1.4109705 -1.4109705 2.5890295 2.5890295
99 100 101 102 103 104 105
-1.1616597 -1.1616597 0.5890295 -0.1616597 -1.1616597 0.8383403 -0.4288542
> postscript(file="/var/wessaorg/rcomp/tmp/642l21323987841.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 105
Frequency = 1
lag(myerror, k = 1) myerror
0 2.5890295 NA
1 -1.4109705 2.5890295
2 2.5711458 -1.4109705
3 -1.1616597 2.5711458
4 -1.4109705 -1.1616597
5 -1.4109705 -1.4109705
6 -1.4109705 -1.4109705
7 -1.4288542 -1.4109705
8 2.5711458 -1.4288542
9 -0.4109705 2.5711458
10 2.8383403 -0.4109705
11 -1.4288542 2.8383403
12 0.5711458 -1.4288542
13 -1.4288542 0.5711458
14 -1.1616597 -1.4288542
15 -1.4109705 -1.1616597
16 -0.4109705 -1.4109705
17 -1.4109705 -0.4109705
18 -1.4288542 -1.4109705
19 0.8383403 -1.4288542
20 0.5890295 0.8383403
21 -0.4109705 0.5890295
22 0.8383403 -0.4109705
23 -1.1437761 0.8383403
24 1.5890295 -1.1437761
25 -1.1437761 1.5890295
26 -1.4109705 -1.1437761
27 -1.1616597 -1.4109705
28 -0.1616597 -1.1616597
29 -1.4109705 -0.1616597
30 -1.1437761 -1.4109705
31 2.5890295 -1.1437761
32 -1.1616597 2.5890295
33 -0.1616597 -1.1616597
34 -1.1616597 -0.1616597
35 -1.4109705 -1.1616597
36 2.5890295 -1.4109705
37 -0.4288542 2.5890295
38 -1.4288542 -0.4288542
39 2.5890295 -1.4288542
40 -1.4109705 2.5890295
41 2.5890295 -1.4109705
42 -1.4109705 2.5890295
43 -1.4109705 -1.4109705
44 -1.1616597 -1.4109705
45 2.5711458 -1.1616597
46 -1.4288542 2.5711458
47 -1.4109705 -1.4288542
48 2.5890295 -1.4109705
49 2.8383403 2.5890295
50 -1.4288542 2.8383403
51 -0.4109705 -1.4288542
52 -1.4288542 -0.4109705
53 2.8383403 -1.4288542
54 -1.4288542 2.8383403
55 0.5711458 -1.4288542
56 -1.4288542 0.5711458
57 2.5711458 -1.4288542
58 2.8383403 2.5711458
59 -1.1616597 2.8383403
60 -1.4288542 -1.1616597
61 2.5890295 -1.4288542
62 -1.4109705 2.5890295
63 -1.1437761 -1.4109705
64 0.8383403 -1.1437761
65 -1.4288542 0.8383403
66 -1.4288542 -1.4288542
67 -1.1616597 -1.4288542
68 2.5890295 -1.1616597
69 2.5890295 2.5890295
70 0.5711458 2.5890295
71 -1.4288542 0.5711458
72 -1.4288542 -1.4288542
73 2.5711458 -1.4288542
74 -1.4109705 2.5711458
75 -1.1437761 -1.4109705
76 -0.1616597 -1.1437761
77 0.5890295 -0.1616597
78 -1.1437761 0.5890295
79 0.5890295 -1.1437761
80 -1.1616597 0.5890295
81 2.8383403 -1.1616597
82 2.8383403 2.8383403
83 -1.1437761 2.8383403
84 -1.1616597 -1.1437761
85 2.8383403 -1.1616597
86 -1.1437761 2.8383403
87 2.5890295 -1.1437761
88 0.8383403 2.5890295
89 0.8383403 0.8383403
90 -1.4109705 0.8383403
91 -1.4109705 -1.4109705
92 2.5890295 -1.4109705
93 -1.4288542 2.5890295
94 -1.4109705 -1.4288542
95 -1.4109705 -1.4109705
96 2.5890295 -1.4109705
97 2.5890295 2.5890295
98 -1.1616597 2.5890295
99 -1.1616597 -1.1616597
100 0.5890295 -1.1616597
101 -0.1616597 0.5890295
102 -1.1616597 -0.1616597
103 0.8383403 -1.1616597
104 -0.4288542 0.8383403
105 NA -0.4288542
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.4109705 2.5890295
[2,] 2.5711458 -1.4109705
[3,] -1.1616597 2.5711458
[4,] -1.4109705 -1.1616597
[5,] -1.4109705 -1.4109705
[6,] -1.4109705 -1.4109705
[7,] -1.4288542 -1.4109705
[8,] 2.5711458 -1.4288542
[9,] -0.4109705 2.5711458
[10,] 2.8383403 -0.4109705
[11,] -1.4288542 2.8383403
[12,] 0.5711458 -1.4288542
[13,] -1.4288542 0.5711458
[14,] -1.1616597 -1.4288542
[15,] -1.4109705 -1.1616597
[16,] -0.4109705 -1.4109705
[17,] -1.4109705 -0.4109705
[18,] -1.4288542 -1.4109705
[19,] 0.8383403 -1.4288542
[20,] 0.5890295 0.8383403
[21,] -0.4109705 0.5890295
[22,] 0.8383403 -0.4109705
[23,] -1.1437761 0.8383403
[24,] 1.5890295 -1.1437761
[25,] -1.1437761 1.5890295
[26,] -1.4109705 -1.1437761
[27,] -1.1616597 -1.4109705
[28,] -0.1616597 -1.1616597
[29,] -1.4109705 -0.1616597
[30,] -1.1437761 -1.4109705
[31,] 2.5890295 -1.1437761
[32,] -1.1616597 2.5890295
[33,] -0.1616597 -1.1616597
[34,] -1.1616597 -0.1616597
[35,] -1.4109705 -1.1616597
[36,] 2.5890295 -1.4109705
[37,] -0.4288542 2.5890295
[38,] -1.4288542 -0.4288542
[39,] 2.5890295 -1.4288542
[40,] -1.4109705 2.5890295
[41,] 2.5890295 -1.4109705
[42,] -1.4109705 2.5890295
[43,] -1.4109705 -1.4109705
[44,] -1.1616597 -1.4109705
[45,] 2.5711458 -1.1616597
[46,] -1.4288542 2.5711458
[47,] -1.4109705 -1.4288542
[48,] 2.5890295 -1.4109705
[49,] 2.8383403 2.5890295
[50,] -1.4288542 2.8383403
[51,] -0.4109705 -1.4288542
[52,] -1.4288542 -0.4109705
[53,] 2.8383403 -1.4288542
[54,] -1.4288542 2.8383403
[55,] 0.5711458 -1.4288542
[56,] -1.4288542 0.5711458
[57,] 2.5711458 -1.4288542
[58,] 2.8383403 2.5711458
[59,] -1.1616597 2.8383403
[60,] -1.4288542 -1.1616597
[61,] 2.5890295 -1.4288542
[62,] -1.4109705 2.5890295
[63,] -1.1437761 -1.4109705
[64,] 0.8383403 -1.1437761
[65,] -1.4288542 0.8383403
[66,] -1.4288542 -1.4288542
[67,] -1.1616597 -1.4288542
[68,] 2.5890295 -1.1616597
[69,] 2.5890295 2.5890295
[70,] 0.5711458 2.5890295
[71,] -1.4288542 0.5711458
[72,] -1.4288542 -1.4288542
[73,] 2.5711458 -1.4288542
[74,] -1.4109705 2.5711458
[75,] -1.1437761 -1.4109705
[76,] -0.1616597 -1.1437761
[77,] 0.5890295 -0.1616597
[78,] -1.1437761 0.5890295
[79,] 0.5890295 -1.1437761
[80,] -1.1616597 0.5890295
[81,] 2.8383403 -1.1616597
[82,] 2.8383403 2.8383403
[83,] -1.1437761 2.8383403
[84,] -1.1616597 -1.1437761
[85,] 2.8383403 -1.1616597
[86,] -1.1437761 2.8383403
[87,] 2.5890295 -1.1437761
[88,] 0.8383403 2.5890295
[89,] 0.8383403 0.8383403
[90,] -1.4109705 0.8383403
[91,] -1.4109705 -1.4109705
[92,] 2.5890295 -1.4109705
[93,] -1.4288542 2.5890295
[94,] -1.4109705 -1.4288542
[95,] -1.4109705 -1.4109705
[96,] 2.5890295 -1.4109705
[97,] 2.5890295 2.5890295
[98,] -1.1616597 2.5890295
[99,] -1.1616597 -1.1616597
[100,] 0.5890295 -1.1616597
[101,] -0.1616597 0.5890295
[102,] -1.1616597 -0.1616597
[103,] 0.8383403 -1.1616597
[104,] -0.4288542 0.8383403
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.4109705 2.5890295
2 2.5711458 -1.4109705
3 -1.1616597 2.5711458
4 -1.4109705 -1.1616597
5 -1.4109705 -1.4109705
6 -1.4109705 -1.4109705
7 -1.4288542 -1.4109705
8 2.5711458 -1.4288542
9 -0.4109705 2.5711458
10 2.8383403 -0.4109705
11 -1.4288542 2.8383403
12 0.5711458 -1.4288542
13 -1.4288542 0.5711458
14 -1.1616597 -1.4288542
15 -1.4109705 -1.1616597
16 -0.4109705 -1.4109705
17 -1.4109705 -0.4109705
18 -1.4288542 -1.4109705
19 0.8383403 -1.4288542
20 0.5890295 0.8383403
21 -0.4109705 0.5890295
22 0.8383403 -0.4109705
23 -1.1437761 0.8383403
24 1.5890295 -1.1437761
25 -1.1437761 1.5890295
26 -1.4109705 -1.1437761
27 -1.1616597 -1.4109705
28 -0.1616597 -1.1616597
29 -1.4109705 -0.1616597
30 -1.1437761 -1.4109705
31 2.5890295 -1.1437761
32 -1.1616597 2.5890295
33 -0.1616597 -1.1616597
34 -1.1616597 -0.1616597
35 -1.4109705 -1.1616597
36 2.5890295 -1.4109705
37 -0.4288542 2.5890295
38 -1.4288542 -0.4288542
39 2.5890295 -1.4288542
40 -1.4109705 2.5890295
41 2.5890295 -1.4109705
42 -1.4109705 2.5890295
43 -1.4109705 -1.4109705
44 -1.1616597 -1.4109705
45 2.5711458 -1.1616597
46 -1.4288542 2.5711458
47 -1.4109705 -1.4288542
48 2.5890295 -1.4109705
49 2.8383403 2.5890295
50 -1.4288542 2.8383403
51 -0.4109705 -1.4288542
52 -1.4288542 -0.4109705
53 2.8383403 -1.4288542
54 -1.4288542 2.8383403
55 0.5711458 -1.4288542
56 -1.4288542 0.5711458
57 2.5711458 -1.4288542
58 2.8383403 2.5711458
59 -1.1616597 2.8383403
60 -1.4288542 -1.1616597
61 2.5890295 -1.4288542
62 -1.4109705 2.5890295
63 -1.1437761 -1.4109705
64 0.8383403 -1.1437761
65 -1.4288542 0.8383403
66 -1.4288542 -1.4288542
67 -1.1616597 -1.4288542
68 2.5890295 -1.1616597
69 2.5890295 2.5890295
70 0.5711458 2.5890295
71 -1.4288542 0.5711458
72 -1.4288542 -1.4288542
73 2.5711458 -1.4288542
74 -1.4109705 2.5711458
75 -1.1437761 -1.4109705
76 -0.1616597 -1.1437761
77 0.5890295 -0.1616597
78 -1.1437761 0.5890295
79 0.5890295 -1.1437761
80 -1.1616597 0.5890295
81 2.8383403 -1.1616597
82 2.8383403 2.8383403
83 -1.1437761 2.8383403
84 -1.1616597 -1.1437761
85 2.8383403 -1.1616597
86 -1.1437761 2.8383403
87 2.5890295 -1.1437761
88 0.8383403 2.5890295
89 0.8383403 0.8383403
90 -1.4109705 0.8383403
91 -1.4109705 -1.4109705
92 2.5890295 -1.4109705
93 -1.4288542 2.5890295
94 -1.4109705 -1.4288542
95 -1.4109705 -1.4109705
96 2.5890295 -1.4109705
97 2.5890295 2.5890295
98 -1.1616597 2.5890295
99 -1.1616597 -1.1616597
100 0.5890295 -1.1616597
101 -0.1616597 0.5890295
102 -1.1616597 -0.1616597
103 0.8383403 -1.1616597
104 -0.4288542 0.8383403
> 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/wessaorg/rcomp/tmp/7wp4p1323987841.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/8404w1323987841.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/9p6gg1323987841.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/10kc0t1323987841.ps",horizontal=F,onefile=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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11gvlw1323987841.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/wessaorg/rcomp/tmp/12i44u1323987841.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/wessaorg/rcomp/tmp/13mpqq1323987841.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/wessaorg/rcomp/tmp/14rbto1323987841.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/wessaorg/rcomp/tmp/151xd51323987841.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/wessaorg/rcomp/tmp/16yi1n1323987841.tab")
+ }
>
> try(system("convert tmp/1ro1k1323987841.ps tmp/1ro1k1323987841.png",intern=TRUE))
character(0)
> try(system("convert tmp/2rewj1323987841.ps tmp/2rewj1323987841.png",intern=TRUE))
character(0)
> try(system("convert tmp/383vs1323987841.ps tmp/383vs1323987841.png",intern=TRUE))
character(0)
> try(system("convert tmp/4gni81323987841.ps tmp/4gni81323987841.png",intern=TRUE))
character(0)
> try(system("convert tmp/5chkv1323987841.ps tmp/5chkv1323987841.png",intern=TRUE))
character(0)
> try(system("convert tmp/642l21323987841.ps tmp/642l21323987841.png",intern=TRUE))
character(0)
> try(system("convert tmp/7wp4p1323987841.ps tmp/7wp4p1323987841.png",intern=TRUE))
character(0)
> try(system("convert tmp/8404w1323987841.ps tmp/8404w1323987841.png",intern=TRUE))
character(0)
> try(system("convert tmp/9p6gg1323987841.ps tmp/9p6gg1323987841.png",intern=TRUE))
character(0)
> try(system("convert tmp/10kc0t1323987841.ps tmp/10kc0t1323987841.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.691 0.603 4.400