R version 2.12.0 (2010-10-15)
Copyright (C) 2010 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(921365,0,987921,0,1132614,0,1332224,0,1418133,0,1411549,0,1695920,0,1636173,0,1539653,0,1395314,0,1127575,0,1036076,0,989236,0,1008380,0,1207763,0,1368839,0,1469798,0,1498721,0,1761769,0,1653214,0,1599104,0,1421179,0,1163995,0,1037735,0,1015407,0,1039210,0,1258049,0,1469445,0,1552346,0,1549144,0,1785895,0,1662335,0,1629440,0,1467430,0,1202209,0,1076982,0,1039367,0,1063449,0,1335135,0,1491602,0,1591972,0,1641248,0,1898849,0,1798580,0,1762444,0,1622044,0,1368955,0,1262973,0,1195650,0,1269530,0,1479279,0,1607819,0,1712466,0,1721766,0,1949843,1,1821326,1,1757802,1,1590367,1,1260647,1,1149235,1,1016367,1,1027885,1,1262159,1,1520854,1,1544144,1,1564709,1,1821776,1,1741365,1,1623386,1,1498658,1,1241822,1,1136029,1),dim=c(2,72),dimnames=list(c('passagiers','dummy'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('passagiers','dummy'),1:72))
> 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
> 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
passagiers dummy M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 921365 0 1 0 0 0 0 0 0 0 0 0 0
2 987921 0 0 1 0 0 0 0 0 0 0 0 0
3 1132614 0 0 0 1 0 0 0 0 0 0 0 0
4 1332224 0 0 0 0 1 0 0 0 0 0 0 0
5 1418133 0 0 0 0 0 1 0 0 0 0 0 0
6 1411549 0 0 0 0 0 0 1 0 0 0 0 0
7 1695920 0 0 0 0 0 0 0 1 0 0 0 0
8 1636173 0 0 0 0 0 0 0 0 1 0 0 0
9 1539653 0 0 0 0 0 0 0 0 0 1 0 0
10 1395314 0 0 0 0 0 0 0 0 0 0 1 0
11 1127575 0 0 0 0 0 0 0 0 0 0 0 1
12 1036076 0 0 0 0 0 0 0 0 0 0 0 0
13 989236 0 1 0 0 0 0 0 0 0 0 0 0
14 1008380 0 0 1 0 0 0 0 0 0 0 0 0
15 1207763 0 0 0 1 0 0 0 0 0 0 0 0
16 1368839 0 0 0 0 1 0 0 0 0 0 0 0
17 1469798 0 0 0 0 0 1 0 0 0 0 0 0
18 1498721 0 0 0 0 0 0 1 0 0 0 0 0
19 1761769 0 0 0 0 0 0 0 1 0 0 0 0
20 1653214 0 0 0 0 0 0 0 0 1 0 0 0
21 1599104 0 0 0 0 0 0 0 0 0 1 0 0
22 1421179 0 0 0 0 0 0 0 0 0 0 1 0
23 1163995 0 0 0 0 0 0 0 0 0 0 0 1
24 1037735 0 0 0 0 0 0 0 0 0 0 0 0
25 1015407 0 1 0 0 0 0 0 0 0 0 0 0
26 1039210 0 0 1 0 0 0 0 0 0 0 0 0
27 1258049 0 0 0 1 0 0 0 0 0 0 0 0
28 1469445 0 0 0 0 1 0 0 0 0 0 0 0
29 1552346 0 0 0 0 0 1 0 0 0 0 0 0
30 1549144 0 0 0 0 0 0 1 0 0 0 0 0
31 1785895 0 0 0 0 0 0 0 1 0 0 0 0
32 1662335 0 0 0 0 0 0 0 0 1 0 0 0
33 1629440 0 0 0 0 0 0 0 0 0 1 0 0
34 1467430 0 0 0 0 0 0 0 0 0 0 1 0
35 1202209 0 0 0 0 0 0 0 0 0 0 0 1
36 1076982 0 0 0 0 0 0 0 0 0 0 0 0
37 1039367 0 1 0 0 0 0 0 0 0 0 0 0
38 1063449 0 0 1 0 0 0 0 0 0 0 0 0
39 1335135 0 0 0 1 0 0 0 0 0 0 0 0
40 1491602 0 0 0 0 1 0 0 0 0 0 0 0
41 1591972 0 0 0 0 0 1 0 0 0 0 0 0
42 1641248 0 0 0 0 0 0 1 0 0 0 0 0
43 1898849 0 0 0 0 0 0 0 1 0 0 0 0
44 1798580 0 0 0 0 0 0 0 0 1 0 0 0
45 1762444 0 0 0 0 0 0 0 0 0 1 0 0
46 1622044 0 0 0 0 0 0 0 0 0 0 1 0
47 1368955 0 0 0 0 0 0 0 0 0 0 0 1
48 1262973 0 0 0 0 0 0 0 0 0 0 0 0
49 1195650 0 1 0 0 0 0 0 0 0 0 0 0
50 1269530 0 0 1 0 0 0 0 0 0 0 0 0
51 1479279 0 0 0 1 0 0 0 0 0 0 0 0
52 1607819 0 0 0 0 1 0 0 0 0 0 0 0
53 1712466 0 0 0 0 0 1 0 0 0 0 0 0
54 1721766 0 0 0 0 0 0 1 0 0 0 0 0
55 1949843 1 0 0 0 0 0 0 1 0 0 0 0
56 1821326 1 0 0 0 0 0 0 0 1 0 0 0
57 1757802 1 0 0 0 0 0 0 0 0 1 0 0
58 1590367 1 0 0 0 0 0 0 0 0 0 1 0
59 1260647 1 0 0 0 0 0 0 0 0 0 0 1
60 1149235 1 0 0 0 0 0 0 0 0 0 0 0
61 1016367 1 1 0 0 0 0 0 0 0 0 0 0
62 1027885 1 0 1 0 0 0 0 0 0 0 0 0
63 1262159 1 0 0 1 0 0 0 0 0 0 0 0
64 1520854 1 0 0 0 1 0 0 0 0 0 0 0
65 1544144 1 0 0 0 0 1 0 0 0 0 0 0
66 1564709 1 0 0 0 0 0 1 0 0 0 0 0
67 1821776 1 0 0 0 0 0 0 1 0 0 0 0
68 1741365 1 0 0 0 0 0 0 0 1 0 0 0
69 1623386 1 0 0 0 0 0 0 0 0 1 0 0
70 1498658 1 0 0 0 0 0 0 0 0 0 1 0
71 1241822 1 0 0 0 0 0 0 0 0 0 0 1
72 1136029 1 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
1103434 39212 -80404 -43907 169197 355161
M5 M6 M7 M8 M9 M10
438173 454553 702504 602327 535466 382660
M11
111029
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-146439 -60667 -19354 53399 210003
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1103434 39892 27.660 < 2e-16 ***
dummy 39212 26432 1.483 0.14327
M1 -80404 55199 -1.457 0.15052
M2 -43907 55199 -0.795 0.42955
M3 169197 55199 3.065 0.00328 **
M4 355161 55199 6.434 2.44e-08 ***
M5 438173 55199 7.938 6.96e-11 ***
M6 454553 55199 8.235 2.20e-11 ***
M7 702504 55023 12.767 < 2e-16 ***
M8 602327 55023 10.947 7.78e-16 ***
M9 535466 55023 9.732 7.05e-14 ***
M10 382660 55023 6.955 3.23e-09 ***
M11 111029 55023 2.018 0.04816 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 95300 on 59 degrees of freedom
Multiple R-squared: 0.897, Adjusted R-squared: 0.876
F-statistic: 42.81 on 12 and 59 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.151578992 0.303157984 0.84842101
[2,] 0.091987536 0.183975072 0.90801246
[3,] 0.094687372 0.189374744 0.90531263
[4,] 0.070988561 0.141977122 0.92901144
[5,] 0.037581858 0.075163716 0.96241814
[6,] 0.027119349 0.054238698 0.97288065
[7,] 0.015875711 0.031751422 0.98412429
[8,] 0.009690662 0.019381323 0.99030934
[9,] 0.005378191 0.010756382 0.99462181
[10,] 0.004747337 0.009494673 0.99525266
[11,] 0.003161290 0.006322580 0.99683871
[12,] 0.005639844 0.011279689 0.99436016
[13,] 0.015302610 0.030605220 0.98469739
[14,] 0.024560288 0.049120576 0.97543971
[15,] 0.031907926 0.063815851 0.96809207
[16,] 0.032514951 0.065029901 0.96748505
[17,] 0.033687027 0.067374053 0.96631297
[18,] 0.036664986 0.073329973 0.96333501
[19,] 0.048131989 0.096263979 0.95186801
[20,] 0.061429506 0.122859011 0.93857049
[21,] 0.095023677 0.190047354 0.90497632
[22,] 0.114381009 0.228762018 0.88561899
[23,] 0.153329980 0.306659960 0.84667002
[24,] 0.266960789 0.533921577 0.73303921
[25,] 0.395259403 0.790518806 0.60474060
[26,] 0.486748920 0.973497839 0.51325108
[27,] 0.592843847 0.814312305 0.40715615
[28,] 0.737897449 0.524205102 0.26210255
[29,] 0.855674143 0.288651714 0.14432586
[30,] 0.899469822 0.201060356 0.10053018
[31,] 0.936009700 0.127980601 0.06399030
[32,] 0.945172698 0.109654605 0.05482730
[33,] 0.951145798 0.097708403 0.04885420
[34,] 0.942199237 0.115601526 0.05780076
[35,] 0.951112237 0.097775526 0.04888776
[36,] 0.951007286 0.097985428 0.04899271
[37,] 0.935575808 0.128848385 0.06442419
[38,] 0.896551869 0.206896263 0.10344813
[39,] 0.831631530 0.336736940 0.16836847
[40,] 0.834757840 0.330484321 0.16524216
[41,] 0.754804520 0.490390960 0.24519548
> postscript(file="/var/www/rcomp/tmp/1n9eo1293456971.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/rcomp/tmp/2f0d91293456971.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/rcomp/tmp/3f0d91293456971.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/rcomp/tmp/4f0d91293456971.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/rcomp/tmp/5q9uu1293456971.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 = 72
Frequency = 1
1 2 3 4 5 6
-101665.075 -71606.241 -140017.241 -126371.241 -123474.908 -146438.575
7 8 9 10 11 12
-110018.150 -69588.650 -99247.983 -90780.816 -86888.316 -67358.483
13 14 15 16 17 18
-33794.075 -51147.241 -64868.241 -89756.241 -71809.908 -59266.575
19 20 21 22 23 24
-44169.150 -52547.650 -39796.983 -64915.816 -50468.316 -65699.483
25 26 27 28 29 30
-7623.075 -20317.241 -14582.241 10849.759 10738.092 -8843.575
31 32 33 34 35 36
-20043.150 -43426.650 -9460.983 -18664.816 -12254.316 -26452.483
37 38 39 40 41 42
16336.925 3921.759 62503.759 33006.759 50364.092 83260.425
43 44 45 46 47 48
92910.850 92818.350 123543.017 135949.184 154491.684 159538.517
49 50 51 52 53 54
172619.925 210002.759 206647.759 149223.759 170858.092 163778.425
55 56 57 58 59 60
104693.299 76352.799 79689.466 65060.632 6972.132 6588.966
61 62 63 64 65 66
-45874.626 -70853.793 -49683.793 23047.207 -36675.459 -32490.126
67 68 69 70 71 72
-23373.701 -3608.201 -54726.534 -26648.368 -11852.868 -6617.034
> postscript(file="/var/www/rcomp/tmp/6q9uu1293456971.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 -101665.075 NA
1 -71606.241 -101665.075
2 -140017.241 -71606.241
3 -126371.241 -140017.241
4 -123474.908 -126371.241
5 -146438.575 -123474.908
6 -110018.150 -146438.575
7 -69588.650 -110018.150
8 -99247.983 -69588.650
9 -90780.816 -99247.983
10 -86888.316 -90780.816
11 -67358.483 -86888.316
12 -33794.075 -67358.483
13 -51147.241 -33794.075
14 -64868.241 -51147.241
15 -89756.241 -64868.241
16 -71809.908 -89756.241
17 -59266.575 -71809.908
18 -44169.150 -59266.575
19 -52547.650 -44169.150
20 -39796.983 -52547.650
21 -64915.816 -39796.983
22 -50468.316 -64915.816
23 -65699.483 -50468.316
24 -7623.075 -65699.483
25 -20317.241 -7623.075
26 -14582.241 -20317.241
27 10849.759 -14582.241
28 10738.092 10849.759
29 -8843.575 10738.092
30 -20043.150 -8843.575
31 -43426.650 -20043.150
32 -9460.983 -43426.650
33 -18664.816 -9460.983
34 -12254.316 -18664.816
35 -26452.483 -12254.316
36 16336.925 -26452.483
37 3921.759 16336.925
38 62503.759 3921.759
39 33006.759 62503.759
40 50364.092 33006.759
41 83260.425 50364.092
42 92910.850 83260.425
43 92818.350 92910.850
44 123543.017 92818.350
45 135949.184 123543.017
46 154491.684 135949.184
47 159538.517 154491.684
48 172619.925 159538.517
49 210002.759 172619.925
50 206647.759 210002.759
51 149223.759 206647.759
52 170858.092 149223.759
53 163778.425 170858.092
54 104693.299 163778.425
55 76352.799 104693.299
56 79689.466 76352.799
57 65060.632 79689.466
58 6972.132 65060.632
59 6588.966 6972.132
60 -45874.626 6588.966
61 -70853.793 -45874.626
62 -49683.793 -70853.793
63 23047.207 -49683.793
64 -36675.459 23047.207
65 -32490.126 -36675.459
66 -23373.701 -32490.126
67 -3608.201 -23373.701
68 -54726.534 -3608.201
69 -26648.368 -54726.534
70 -11852.868 -26648.368
71 -6617.034 -11852.868
72 NA -6617.034
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -71606.241 -101665.075
[2,] -140017.241 -71606.241
[3,] -126371.241 -140017.241
[4,] -123474.908 -126371.241
[5,] -146438.575 -123474.908
[6,] -110018.150 -146438.575
[7,] -69588.650 -110018.150
[8,] -99247.983 -69588.650
[9,] -90780.816 -99247.983
[10,] -86888.316 -90780.816
[11,] -67358.483 -86888.316
[12,] -33794.075 -67358.483
[13,] -51147.241 -33794.075
[14,] -64868.241 -51147.241
[15,] -89756.241 -64868.241
[16,] -71809.908 -89756.241
[17,] -59266.575 -71809.908
[18,] -44169.150 -59266.575
[19,] -52547.650 -44169.150
[20,] -39796.983 -52547.650
[21,] -64915.816 -39796.983
[22,] -50468.316 -64915.816
[23,] -65699.483 -50468.316
[24,] -7623.075 -65699.483
[25,] -20317.241 -7623.075
[26,] -14582.241 -20317.241
[27,] 10849.759 -14582.241
[28,] 10738.092 10849.759
[29,] -8843.575 10738.092
[30,] -20043.150 -8843.575
[31,] -43426.650 -20043.150
[32,] -9460.983 -43426.650
[33,] -18664.816 -9460.983
[34,] -12254.316 -18664.816
[35,] -26452.483 -12254.316
[36,] 16336.925 -26452.483
[37,] 3921.759 16336.925
[38,] 62503.759 3921.759
[39,] 33006.759 62503.759
[40,] 50364.092 33006.759
[41,] 83260.425 50364.092
[42,] 92910.850 83260.425
[43,] 92818.350 92910.850
[44,] 123543.017 92818.350
[45,] 135949.184 123543.017
[46,] 154491.684 135949.184
[47,] 159538.517 154491.684
[48,] 172619.925 159538.517
[49,] 210002.759 172619.925
[50,] 206647.759 210002.759
[51,] 149223.759 206647.759
[52,] 170858.092 149223.759
[53,] 163778.425 170858.092
[54,] 104693.299 163778.425
[55,] 76352.799 104693.299
[56,] 79689.466 76352.799
[57,] 65060.632 79689.466
[58,] 6972.132 65060.632
[59,] 6588.966 6972.132
[60,] -45874.626 6588.966
[61,] -70853.793 -45874.626
[62,] -49683.793 -70853.793
[63,] 23047.207 -49683.793
[64,] -36675.459 23047.207
[65,] -32490.126 -36675.459
[66,] -23373.701 -32490.126
[67,] -3608.201 -23373.701
[68,] -54726.534 -3608.201
[69,] -26648.368 -54726.534
[70,] -11852.868 -26648.368
[71,] -6617.034 -11852.868
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -71606.241 -101665.075
2 -140017.241 -71606.241
3 -126371.241 -140017.241
4 -123474.908 -126371.241
5 -146438.575 -123474.908
6 -110018.150 -146438.575
7 -69588.650 -110018.150
8 -99247.983 -69588.650
9 -90780.816 -99247.983
10 -86888.316 -90780.816
11 -67358.483 -86888.316
12 -33794.075 -67358.483
13 -51147.241 -33794.075
14 -64868.241 -51147.241
15 -89756.241 -64868.241
16 -71809.908 -89756.241
17 -59266.575 -71809.908
18 -44169.150 -59266.575
19 -52547.650 -44169.150
20 -39796.983 -52547.650
21 -64915.816 -39796.983
22 -50468.316 -64915.816
23 -65699.483 -50468.316
24 -7623.075 -65699.483
25 -20317.241 -7623.075
26 -14582.241 -20317.241
27 10849.759 -14582.241
28 10738.092 10849.759
29 -8843.575 10738.092
30 -20043.150 -8843.575
31 -43426.650 -20043.150
32 -9460.983 -43426.650
33 -18664.816 -9460.983
34 -12254.316 -18664.816
35 -26452.483 -12254.316
36 16336.925 -26452.483
37 3921.759 16336.925
38 62503.759 3921.759
39 33006.759 62503.759
40 50364.092 33006.759
41 83260.425 50364.092
42 92910.850 83260.425
43 92818.350 92910.850
44 123543.017 92818.350
45 135949.184 123543.017
46 154491.684 135949.184
47 159538.517 154491.684
48 172619.925 159538.517
49 210002.759 172619.925
50 206647.759 210002.759
51 149223.759 206647.759
52 170858.092 149223.759
53 163778.425 170858.092
54 104693.299 163778.425
55 76352.799 104693.299
56 79689.466 76352.799
57 65060.632 79689.466
58 6972.132 65060.632
59 6588.966 6972.132
60 -45874.626 6588.966
61 -70853.793 -45874.626
62 -49683.793 -70853.793
63 23047.207 -49683.793
64 -36675.459 23047.207
65 -32490.126 -36675.459
66 -23373.701 -32490.126
67 -3608.201 -23373.701
68 -54726.534 -3608.201
69 -26648.368 -54726.534
70 -11852.868 -26648.368
71 -6617.034 -11852.868
> 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/rcomp/tmp/7j1uf1293456971.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/rcomp/tmp/8j1uf1293456971.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/rcomp/tmp/9tsaz1293456971.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/rcomp/tmp/10tsaz1293456971.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/11xsrn1293456971.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/rcomp/tmp/12ibqt1293456971.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/rcomp/tmp/13hm7i1293456972.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/rcomp/tmp/14sdo31293456972.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/rcomp/tmp/15dwn91293456972.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/rcomp/tmp/1695lh1293456972.tab")
+ }
>
> try(system("convert tmp/1n9eo1293456971.ps tmp/1n9eo1293456971.png",intern=TRUE))
character(0)
> try(system("convert tmp/2f0d91293456971.ps tmp/2f0d91293456971.png",intern=TRUE))
character(0)
> try(system("convert tmp/3f0d91293456971.ps tmp/3f0d91293456971.png",intern=TRUE))
character(0)
> try(system("convert tmp/4f0d91293456971.ps tmp/4f0d91293456971.png",intern=TRUE))
character(0)
> try(system("convert tmp/5q9uu1293456971.ps tmp/5q9uu1293456971.png",intern=TRUE))
character(0)
> try(system("convert tmp/6q9uu1293456971.ps tmp/6q9uu1293456971.png",intern=TRUE))
character(0)
> try(system("convert tmp/7j1uf1293456971.ps tmp/7j1uf1293456971.png",intern=TRUE))
character(0)
> try(system("convert tmp/8j1uf1293456971.ps tmp/8j1uf1293456971.png",intern=TRUE))
character(0)
> try(system("convert tmp/9tsaz1293456971.ps tmp/9tsaz1293456971.png",intern=TRUE))
character(0)
> try(system("convert tmp/10tsaz1293456971.ps tmp/10tsaz1293456971.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.330 1.640 4.956