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(4940,1,3924,1,3927,1,4535,1,3446,1,3016,1,4934,1,2743,1,3242,1,6662,1,3262,1,3381,1,7144,1,3803,1,3684,1,6759,1,3386,1,3066,1,5538,1,2940,1,3215,1,7023,1,3443,1,3712,1,7475,1,4137,1,3491,1,7019,1,3908,1,3402,1,5604,1,3222,1,3636,1,7123,1,4368,1,4092,1,8377,1,4595,1,4188,1,6988,1,4218,1,3655,1,6211,1,3622,1,3841,1,8510,1,4627,1,4582,1,8967,1,4928,1,4809,1,7917,1,4790,1,4065,1,7290,1,4670,1,3561,1,5149,1,6880,1,6981,1,8454,1,4960,1,4670,1,7638,1,4560,1,3980,0,6825,0,3939,0,4079,0,8117,0,5121,0,5167,0,7960,0,4670,0,4397,0,7191,0,4293,0,3747,0,6425,0,3709,0,3840,0,7642,0,4821,0,4865,0),dim=c(2,84),dimnames=list(c('ASO','Dummy'),1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('ASO','Dummy'),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 = '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
ASO Dummy M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 4940 1 1 0 0 0 0 0 0 0 0 0 0
2 3924 1 0 1 0 0 0 0 0 0 0 0 0
3 3927 1 0 0 1 0 0 0 0 0 0 0 0
4 4535 1 0 0 0 1 0 0 0 0 0 0 0
5 3446 1 0 0 0 0 1 0 0 0 0 0 0
6 3016 1 0 0 0 0 0 1 0 0 0 0 0
7 4934 1 0 0 0 0 0 0 1 0 0 0 0
8 2743 1 0 0 0 0 0 0 0 1 0 0 0
9 3242 1 0 0 0 0 0 0 0 0 1 0 0
10 6662 1 0 0 0 0 0 0 0 0 0 1 0
11 3262 1 0 0 0 0 0 0 0 0 0 0 1
12 3381 1 0 0 0 0 0 0 0 0 0 0 0
13 7144 1 1 0 0 0 0 0 0 0 0 0 0
14 3803 1 0 1 0 0 0 0 0 0 0 0 0
15 3684 1 0 0 1 0 0 0 0 0 0 0 0
16 6759 1 0 0 0 1 0 0 0 0 0 0 0
17 3386 1 0 0 0 0 1 0 0 0 0 0 0
18 3066 1 0 0 0 0 0 1 0 0 0 0 0
19 5538 1 0 0 0 0 0 0 1 0 0 0 0
20 2940 1 0 0 0 0 0 0 0 1 0 0 0
21 3215 1 0 0 0 0 0 0 0 0 1 0 0
22 7023 1 0 0 0 0 0 0 0 0 0 1 0
23 3443 1 0 0 0 0 0 0 0 0 0 0 1
24 3712 1 0 0 0 0 0 0 0 0 0 0 0
25 7475 1 1 0 0 0 0 0 0 0 0 0 0
26 4137 1 0 1 0 0 0 0 0 0 0 0 0
27 3491 1 0 0 1 0 0 0 0 0 0 0 0
28 7019 1 0 0 0 1 0 0 0 0 0 0 0
29 3908 1 0 0 0 0 1 0 0 0 0 0 0
30 3402 1 0 0 0 0 0 1 0 0 0 0 0
31 5604 1 0 0 0 0 0 0 1 0 0 0 0
32 3222 1 0 0 0 0 0 0 0 1 0 0 0
33 3636 1 0 0 0 0 0 0 0 0 1 0 0
34 7123 1 0 0 0 0 0 0 0 0 0 1 0
35 4368 1 0 0 0 0 0 0 0 0 0 0 1
36 4092 1 0 0 0 0 0 0 0 0 0 0 0
37 8377 1 1 0 0 0 0 0 0 0 0 0 0
38 4595 1 0 1 0 0 0 0 0 0 0 0 0
39 4188 1 0 0 1 0 0 0 0 0 0 0 0
40 6988 1 0 0 0 1 0 0 0 0 0 0 0
41 4218 1 0 0 0 0 1 0 0 0 0 0 0
42 3655 1 0 0 0 0 0 1 0 0 0 0 0
43 6211 1 0 0 0 0 0 0 1 0 0 0 0
44 3622 1 0 0 0 0 0 0 0 1 0 0 0
45 3841 1 0 0 0 0 0 0 0 0 1 0 0
46 8510 1 0 0 0 0 0 0 0 0 0 1 0
47 4627 1 0 0 0 0 0 0 0 0 0 0 1
48 4582 1 0 0 0 0 0 0 0 0 0 0 0
49 8967 1 1 0 0 0 0 0 0 0 0 0 0
50 4928 1 0 1 0 0 0 0 0 0 0 0 0
51 4809 1 0 0 1 0 0 0 0 0 0 0 0
52 7917 1 0 0 0 1 0 0 0 0 0 0 0
53 4790 1 0 0 0 0 1 0 0 0 0 0 0
54 4065 1 0 0 0 0 0 1 0 0 0 0 0
55 7290 1 0 0 0 0 0 0 1 0 0 0 0
56 4670 1 0 0 0 0 0 0 0 1 0 0 0
57 3561 1 0 0 0 0 0 0 0 0 1 0 0
58 5149 1 0 0 0 0 0 0 0 0 0 1 0
59 6880 1 0 0 0 0 0 0 0 0 0 0 1
60 6981 1 0 0 0 0 0 0 0 0 0 0 0
61 8454 1 1 0 0 0 0 0 0 0 0 0 0
62 4960 1 0 1 0 0 0 0 0 0 0 0 0
63 4670 1 0 0 1 0 0 0 0 0 0 0 0
64 7638 1 0 0 0 1 0 0 0 0 0 0 0
65 4560 1 0 0 0 0 1 0 0 0 0 0 0
66 3980 0 0 0 0 0 0 1 0 0 0 0 0
67 6825 0 0 0 0 0 0 0 1 0 0 0 0
68 3939 0 0 0 0 0 0 0 0 1 0 0 0
69 4079 0 0 0 0 0 0 0 0 0 1 0 0
70 8117 0 0 0 0 0 0 0 0 0 0 1 0
71 5121 0 0 0 0 0 0 0 0 0 0 0 1
72 5167 0 0 0 0 0 0 0 0 0 0 0 0
73 7960 0 1 0 0 0 0 0 0 0 0 0 0
74 4670 0 0 1 0 0 0 0 0 0 0 0 0
75 4397 0 0 0 1 0 0 0 0 0 0 0 0
76 7191 0 0 0 0 1 0 0 0 0 0 0 0
77 4293 0 0 0 0 0 1 0 0 0 0 0 0
78 3747 0 0 0 0 0 0 1 0 0 0 0 0
79 6425 0 0 0 0 0 0 0 1 0 0 0 0
80 3709 0 0 0 0 0 0 0 0 1 0 0 0
81 3840 0 0 0 0 0 0 0 0 0 1 0 0
82 7642 0 0 0 0 0 0 0 0 0 0 1 0
83 4821 0 0 0 0 0 0 0 0 0 0 0 1
84 4865 0 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Dummy M1 M2 M3 M4
5027.70 -482.78 3002.83 -182.89 -447.32 2249.97
M5 M6 M7 M8 M9 M10
-528.03 -1121.29 1435.29 -1133.57 -1052.29 2492.29
M11
-36.86
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2607.75 -375.46 -36.94 231.77 2436.08
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5027.70 361.06 13.925 < 2e-16 ***
Dummy -482.78 226.06 -2.136 0.03616 *
M1 3002.83 457.85 6.558 7.50e-09 ***
M2 -182.89 457.85 -0.399 0.69076
M3 -447.32 457.85 -0.977 0.33189
M4 2249.97 457.85 4.914 5.53e-06 ***
M5 -528.03 457.85 -1.153 0.25266
M6 -1121.29 456.71 -2.455 0.01654 *
M7 1435.29 456.71 3.143 0.00244 **
M8 -1133.57 456.71 -2.482 0.01543 *
M9 -1052.29 456.71 -2.304 0.02415 *
M10 2492.29 456.71 5.457 6.70e-07 ***
M11 -36.86 456.71 -0.081 0.93591
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 854.4 on 71 degrees of freedom
Multiple R-squared: 0.7707, Adjusted R-squared: 0.7319
F-statistic: 19.88 on 12 and 71 DF, p-value: < 2.2e-16
> 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.9661762 6.764760e-02 3.382380e-02
[2,] 0.9337688 1.324624e-01 6.623121e-02
[3,] 0.8832025 2.335950e-01 1.167975e-01
[4,] 0.8377195 3.245610e-01 1.622805e-01
[5,] 0.7697583 4.604835e-01 2.302417e-01
[6,] 0.6825366 6.349268e-01 3.174634e-01
[7,] 0.5926323 8.147354e-01 4.073677e-01
[8,] 0.5648759 8.702481e-01 4.351241e-01
[9,] 0.5306073 9.387853e-01 4.693927e-01
[10,] 0.6570151 6.859698e-01 3.429849e-01
[11,] 0.5890799 8.218402e-01 4.109201e-01
[12,] 0.5442373 9.115253e-01 4.557627e-01
[13,] 0.6196497 7.607005e-01 3.803503e-01
[14,] 0.5664649 8.670702e-01 4.335351e-01
[15,] 0.5011215 9.977569e-01 4.988785e-01
[16,] 0.4759698 9.519396e-01 5.240302e-01
[17,] 0.4318380 8.636760e-01 5.681620e-01
[18,] 0.3687679 7.375358e-01 6.312321e-01
[19,] 0.3045139 6.090277e-01 6.954861e-01
[20,] 0.3466436 6.932872e-01 6.533564e-01
[21,] 0.3791348 7.582696e-01 6.208652e-01
[22,] 0.5522007 8.955987e-01 4.477993e-01
[23,] 0.5075393 9.849213e-01 4.924607e-01
[24,] 0.4590430 9.180860e-01 5.409570e-01
[25,] 0.4534402 9.068804e-01 5.465598e-01
[26,] 0.4112031 8.224061e-01 5.887969e-01
[27,] 0.3617622 7.235243e-01 6.382378e-01
[28,] 0.3621823 7.243646e-01 6.378177e-01
[29,] 0.3375360 6.750720e-01 6.624640e-01
[30,] 0.2838522 5.677044e-01 7.161478e-01
[31,] 0.4265877 8.531754e-01 5.734123e-01
[32,] 0.4753809 9.507619e-01 5.246191e-01
[33,] 0.5496807 9.006387e-01 4.503193e-01
[34,] 0.6559633 6.880735e-01 3.440367e-01
[35,] 0.6061770 7.876459e-01 3.938230e-01
[36,] 0.5644526 8.710949e-01 4.355474e-01
[37,] 0.5773554 8.452892e-01 4.226446e-01
[38,] 0.5338297 9.323406e-01 4.661703e-01
[39,] 0.4718104 9.436207e-01 5.281896e-01
[40,] 0.4787976 9.575952e-01 5.212024e-01
[41,] 0.4727545 9.455090e-01 5.272455e-01
[42,] 0.4170259 8.340518e-01 5.829741e-01
[43,] 0.9939791 1.204181e-02 6.020906e-03
[44,] 0.9991400 1.719975e-03 8.599874e-04
[45,] 0.9999995 1.013444e-06 5.067222e-07
[46,] 0.9999978 4.383323e-06 2.191661e-06
[47,] 0.9999887 2.266051e-05 1.133026e-05
[48,] 0.9999451 1.097168e-04 5.485840e-05
[49,] 0.9997679 4.642176e-04 2.321088e-04
[50,] 0.9989589 2.082136e-03 1.041068e-03
[51,] 0.9963302 7.339687e-03 3.669843e-03
[52,] 0.9911700 1.765996e-02 8.829981e-03
[53,] 0.9687768 6.244639e-02 3.122320e-02
> postscript(file="/var/www/html/rcomp/tmp/1kn2b1292941575.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/www/html/rcomp/tmp/2kn2b1292941575.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/www/html/rcomp/tmp/3de1w1292941575.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/www/html/rcomp/tmp/4de1w1292941575.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/www/html/rcomp/tmp/5de1w1292941575.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 = 84
Frequency = 1
1 2 3 4 5 6
-2607.74571 -438.03143 -170.60286 -2259.88857 -570.88857 -407.63429
7 8 9 10 11 12
-1046.20571 -668.34857 -250.63429 -375.20571 -1246.06286 -1163.92000
13 14 15 16 17 18
-403.74571 -559.03143 -413.60286 -35.88857 -630.88857 -357.63429
19 20 21 22 23 24
-442.20571 -471.34857 -277.63429 -14.20571 -1065.06286 -832.92000
25 26 27 28 29 30
-72.74571 -225.03143 -606.60286 224.11143 -108.88857 -21.63429
31 32 33 34 35 36
-376.20571 -189.34857 143.36571 85.79429 -140.06286 -452.92000
37 38 39 40 41 42
829.25429 232.96857 90.39714 193.11143 201.11143 231.36571
43 44 45 46 47 48
230.79429 210.65143 348.36571 1472.79429 118.93714 37.08000
49 50 51 52 53 54
1419.25429 565.96857 711.39714 1122.11143 773.11143 641.36571
55 56 57 58 59 60
1309.79429 1258.65143 68.36571 -1888.20571 2371.93714 2436.08000
61 62 63 64 65 66
906.25429 597.96857 572.39714 843.11143 543.11143 73.58571
67 68 69 70 71 72
362.01429 44.87143 103.58571 597.01429 130.15714 139.30000
73 74 75 76 77 78
-70.52571 -174.81143 -183.38286 -86.66857 -206.66857 -159.41429
79 80 81 82 83 84
-37.98571 -185.12857 -135.41429 122.01429 -169.84286 -162.70000
> postscript(file="/var/www/html/rcomp/tmp/6o5iz1292941575.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 = 84
Frequency = 1
lag(myerror, k = 1) myerror
0 -2607.74571 NA
1 -438.03143 -2607.74571
2 -170.60286 -438.03143
3 -2259.88857 -170.60286
4 -570.88857 -2259.88857
5 -407.63429 -570.88857
6 -1046.20571 -407.63429
7 -668.34857 -1046.20571
8 -250.63429 -668.34857
9 -375.20571 -250.63429
10 -1246.06286 -375.20571
11 -1163.92000 -1246.06286
12 -403.74571 -1163.92000
13 -559.03143 -403.74571
14 -413.60286 -559.03143
15 -35.88857 -413.60286
16 -630.88857 -35.88857
17 -357.63429 -630.88857
18 -442.20571 -357.63429
19 -471.34857 -442.20571
20 -277.63429 -471.34857
21 -14.20571 -277.63429
22 -1065.06286 -14.20571
23 -832.92000 -1065.06286
24 -72.74571 -832.92000
25 -225.03143 -72.74571
26 -606.60286 -225.03143
27 224.11143 -606.60286
28 -108.88857 224.11143
29 -21.63429 -108.88857
30 -376.20571 -21.63429
31 -189.34857 -376.20571
32 143.36571 -189.34857
33 85.79429 143.36571
34 -140.06286 85.79429
35 -452.92000 -140.06286
36 829.25429 -452.92000
37 232.96857 829.25429
38 90.39714 232.96857
39 193.11143 90.39714
40 201.11143 193.11143
41 231.36571 201.11143
42 230.79429 231.36571
43 210.65143 230.79429
44 348.36571 210.65143
45 1472.79429 348.36571
46 118.93714 1472.79429
47 37.08000 118.93714
48 1419.25429 37.08000
49 565.96857 1419.25429
50 711.39714 565.96857
51 1122.11143 711.39714
52 773.11143 1122.11143
53 641.36571 773.11143
54 1309.79429 641.36571
55 1258.65143 1309.79429
56 68.36571 1258.65143
57 -1888.20571 68.36571
58 2371.93714 -1888.20571
59 2436.08000 2371.93714
60 906.25429 2436.08000
61 597.96857 906.25429
62 572.39714 597.96857
63 843.11143 572.39714
64 543.11143 843.11143
65 73.58571 543.11143
66 362.01429 73.58571
67 44.87143 362.01429
68 103.58571 44.87143
69 597.01429 103.58571
70 130.15714 597.01429
71 139.30000 130.15714
72 -70.52571 139.30000
73 -174.81143 -70.52571
74 -183.38286 -174.81143
75 -86.66857 -183.38286
76 -206.66857 -86.66857
77 -159.41429 -206.66857
78 -37.98571 -159.41429
79 -185.12857 -37.98571
80 -135.41429 -185.12857
81 122.01429 -135.41429
82 -169.84286 122.01429
83 -162.70000 -169.84286
84 NA -162.70000
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -438.03143 -2607.74571
[2,] -170.60286 -438.03143
[3,] -2259.88857 -170.60286
[4,] -570.88857 -2259.88857
[5,] -407.63429 -570.88857
[6,] -1046.20571 -407.63429
[7,] -668.34857 -1046.20571
[8,] -250.63429 -668.34857
[9,] -375.20571 -250.63429
[10,] -1246.06286 -375.20571
[11,] -1163.92000 -1246.06286
[12,] -403.74571 -1163.92000
[13,] -559.03143 -403.74571
[14,] -413.60286 -559.03143
[15,] -35.88857 -413.60286
[16,] -630.88857 -35.88857
[17,] -357.63429 -630.88857
[18,] -442.20571 -357.63429
[19,] -471.34857 -442.20571
[20,] -277.63429 -471.34857
[21,] -14.20571 -277.63429
[22,] -1065.06286 -14.20571
[23,] -832.92000 -1065.06286
[24,] -72.74571 -832.92000
[25,] -225.03143 -72.74571
[26,] -606.60286 -225.03143
[27,] 224.11143 -606.60286
[28,] -108.88857 224.11143
[29,] -21.63429 -108.88857
[30,] -376.20571 -21.63429
[31,] -189.34857 -376.20571
[32,] 143.36571 -189.34857
[33,] 85.79429 143.36571
[34,] -140.06286 85.79429
[35,] -452.92000 -140.06286
[36,] 829.25429 -452.92000
[37,] 232.96857 829.25429
[38,] 90.39714 232.96857
[39,] 193.11143 90.39714
[40,] 201.11143 193.11143
[41,] 231.36571 201.11143
[42,] 230.79429 231.36571
[43,] 210.65143 230.79429
[44,] 348.36571 210.65143
[45,] 1472.79429 348.36571
[46,] 118.93714 1472.79429
[47,] 37.08000 118.93714
[48,] 1419.25429 37.08000
[49,] 565.96857 1419.25429
[50,] 711.39714 565.96857
[51,] 1122.11143 711.39714
[52,] 773.11143 1122.11143
[53,] 641.36571 773.11143
[54,] 1309.79429 641.36571
[55,] 1258.65143 1309.79429
[56,] 68.36571 1258.65143
[57,] -1888.20571 68.36571
[58,] 2371.93714 -1888.20571
[59,] 2436.08000 2371.93714
[60,] 906.25429 2436.08000
[61,] 597.96857 906.25429
[62,] 572.39714 597.96857
[63,] 843.11143 572.39714
[64,] 543.11143 843.11143
[65,] 73.58571 543.11143
[66,] 362.01429 73.58571
[67,] 44.87143 362.01429
[68,] 103.58571 44.87143
[69,] 597.01429 103.58571
[70,] 130.15714 597.01429
[71,] 139.30000 130.15714
[72,] -70.52571 139.30000
[73,] -174.81143 -70.52571
[74,] -183.38286 -174.81143
[75,] -86.66857 -183.38286
[76,] -206.66857 -86.66857
[77,] -159.41429 -206.66857
[78,] -37.98571 -159.41429
[79,] -185.12857 -37.98571
[80,] -135.41429 -185.12857
[81,] 122.01429 -135.41429
[82,] -169.84286 122.01429
[83,] -162.70000 -169.84286
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -438.03143 -2607.74571
2 -170.60286 -438.03143
3 -2259.88857 -170.60286
4 -570.88857 -2259.88857
5 -407.63429 -570.88857
6 -1046.20571 -407.63429
7 -668.34857 -1046.20571
8 -250.63429 -668.34857
9 -375.20571 -250.63429
10 -1246.06286 -375.20571
11 -1163.92000 -1246.06286
12 -403.74571 -1163.92000
13 -559.03143 -403.74571
14 -413.60286 -559.03143
15 -35.88857 -413.60286
16 -630.88857 -35.88857
17 -357.63429 -630.88857
18 -442.20571 -357.63429
19 -471.34857 -442.20571
20 -277.63429 -471.34857
21 -14.20571 -277.63429
22 -1065.06286 -14.20571
23 -832.92000 -1065.06286
24 -72.74571 -832.92000
25 -225.03143 -72.74571
26 -606.60286 -225.03143
27 224.11143 -606.60286
28 -108.88857 224.11143
29 -21.63429 -108.88857
30 -376.20571 -21.63429
31 -189.34857 -376.20571
32 143.36571 -189.34857
33 85.79429 143.36571
34 -140.06286 85.79429
35 -452.92000 -140.06286
36 829.25429 -452.92000
37 232.96857 829.25429
38 90.39714 232.96857
39 193.11143 90.39714
40 201.11143 193.11143
41 231.36571 201.11143
42 230.79429 231.36571
43 210.65143 230.79429
44 348.36571 210.65143
45 1472.79429 348.36571
46 118.93714 1472.79429
47 37.08000 118.93714
48 1419.25429 37.08000
49 565.96857 1419.25429
50 711.39714 565.96857
51 1122.11143 711.39714
52 773.11143 1122.11143
53 641.36571 773.11143
54 1309.79429 641.36571
55 1258.65143 1309.79429
56 68.36571 1258.65143
57 -1888.20571 68.36571
58 2371.93714 -1888.20571
59 2436.08000 2371.93714
60 906.25429 2436.08000
61 597.96857 906.25429
62 572.39714 597.96857
63 843.11143 572.39714
64 543.11143 843.11143
65 73.58571 543.11143
66 362.01429 73.58571
67 44.87143 362.01429
68 103.58571 44.87143
69 597.01429 103.58571
70 130.15714 597.01429
71 139.30000 130.15714
72 -70.52571 139.30000
73 -174.81143 -70.52571
74 -183.38286 -174.81143
75 -86.66857 -183.38286
76 -206.66857 -86.66857
77 -159.41429 -206.66857
78 -37.98571 -159.41429
79 -185.12857 -37.98571
80 -135.41429 -185.12857
81 122.01429 -135.41429
82 -169.84286 122.01429
83 -162.70000 -169.84286
> 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/7yfik1292941575.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/www/html/rcomp/tmp/8yfik1292941575.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/www/html/rcomp/tmp/9yfik1292941575.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/www/html/rcomp/tmp/10r6z51292941575.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/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/11nyin1292941576.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/12qhgb1292941576.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/13n9wk1292941576.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/148rv81292941576.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/15tste1292941576.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/16fask1292941576.tab")
+ }
>
> try(system("convert tmp/1kn2b1292941575.ps tmp/1kn2b1292941575.png",intern=TRUE))
character(0)
> try(system("convert tmp/2kn2b1292941575.ps tmp/2kn2b1292941575.png",intern=TRUE))
character(0)
> try(system("convert tmp/3de1w1292941575.ps tmp/3de1w1292941575.png",intern=TRUE))
character(0)
> try(system("convert tmp/4de1w1292941575.ps tmp/4de1w1292941575.png",intern=TRUE))
character(0)
> try(system("convert tmp/5de1w1292941575.ps tmp/5de1w1292941575.png",intern=TRUE))
character(0)
> try(system("convert tmp/6o5iz1292941575.ps tmp/6o5iz1292941575.png",intern=TRUE))
character(0)
> try(system("convert tmp/7yfik1292941575.ps tmp/7yfik1292941575.png",intern=TRUE))
character(0)
> try(system("convert tmp/8yfik1292941575.ps tmp/8yfik1292941575.png",intern=TRUE))
character(0)
> try(system("convert tmp/9yfik1292941575.ps tmp/9yfik1292941575.png",intern=TRUE))
character(0)
> try(system("convert tmp/10r6z51292941575.ps tmp/10r6z51292941575.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.789 1.790 9.369