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(921365,987921,1132614,1332224,1418133,1411549,1695920,1636173,1539653,1395314,1127575,1036076,989236,1008380,1207763,1368839,1469798,1498721,1761761,1653214,1599104,1421179,1163995,1037735,1015407,1039210,1258049,1469445,1552346,1549144,1785895,1662335,1629440,1467430,1202209,1076982,1039367,1063449,1335135,1491602,1591972,1641248,1898849,1798580,1762444,1622044,1368955,1262973,1269530,1479279,1607819,1721466,1721766,1949843,1821326,1757802,1590367,1260647,1149235,1016367,1027885,1262159,1520854,1544144,1564709,1821776,1741365,1623386,1498658,1241822,1136029,1035030,1078521,1279431,1171023,1573377,1589514,1859878,1783191,1689849,1619868,1323443,1177481),dim=c(1,83),dimnames=list(c('Passagiers_per_maand'),1:83))
> y <- array(NA,dim=c(1,83),dimnames=list(c('Passagiers_per_maand'),1:83))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
> 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_per_maand M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 921365 1 0 0 0 0 0 0 0 0 0 0 1
2 987921 0 1 0 0 0 0 0 0 0 0 0 2
3 1132614 0 0 1 0 0 0 0 0 0 0 0 3
4 1332224 0 0 0 1 0 0 0 0 0 0 0 4
5 1418133 0 0 0 0 1 0 0 0 0 0 0 5
6 1411549 0 0 0 0 0 1 0 0 0 0 0 6
7 1695920 0 0 0 0 0 0 1 0 0 0 0 7
8 1636173 0 0 0 0 0 0 0 1 0 0 0 8
9 1539653 0 0 0 0 0 0 0 0 1 0 0 9
10 1395314 0 0 0 0 0 0 0 0 0 1 0 10
11 1127575 0 0 0 0 0 0 0 0 0 0 1 11
12 1036076 0 0 0 0 0 0 0 0 0 0 0 12
13 989236 1 0 0 0 0 0 0 0 0 0 0 13
14 1008380 0 1 0 0 0 0 0 0 0 0 0 14
15 1207763 0 0 1 0 0 0 0 0 0 0 0 15
16 1368839 0 0 0 1 0 0 0 0 0 0 0 16
17 1469798 0 0 0 0 1 0 0 0 0 0 0 17
18 1498721 0 0 0 0 0 1 0 0 0 0 0 18
19 1761761 0 0 0 0 0 0 1 0 0 0 0 19
20 1653214 0 0 0 0 0 0 0 1 0 0 0 20
21 1599104 0 0 0 0 0 0 0 0 1 0 0 21
22 1421179 0 0 0 0 0 0 0 0 0 1 0 22
23 1163995 0 0 0 0 0 0 0 0 0 0 1 23
24 1037735 0 0 0 0 0 0 0 0 0 0 0 24
25 1015407 1 0 0 0 0 0 0 0 0 0 0 25
26 1039210 0 1 0 0 0 0 0 0 0 0 0 26
27 1258049 0 0 1 0 0 0 0 0 0 0 0 27
28 1469445 0 0 0 1 0 0 0 0 0 0 0 28
29 1552346 0 0 0 0 1 0 0 0 0 0 0 29
30 1549144 0 0 0 0 0 1 0 0 0 0 0 30
31 1785895 0 0 0 0 0 0 1 0 0 0 0 31
32 1662335 0 0 0 0 0 0 0 1 0 0 0 32
33 1629440 0 0 0 0 0 0 0 0 1 0 0 33
34 1467430 0 0 0 0 0 0 0 0 0 1 0 34
35 1202209 0 0 0 0 0 0 0 0 0 0 1 35
36 1076982 0 0 0 0 0 0 0 0 0 0 0 36
37 1039367 1 0 0 0 0 0 0 0 0 0 0 37
38 1063449 0 1 0 0 0 0 0 0 0 0 0 38
39 1335135 0 0 1 0 0 0 0 0 0 0 0 39
40 1491602 0 0 0 1 0 0 0 0 0 0 0 40
41 1591972 0 0 0 0 1 0 0 0 0 0 0 41
42 1641248 0 0 0 0 0 1 0 0 0 0 0 42
43 1898849 0 0 0 0 0 0 1 0 0 0 0 43
44 1798580 0 0 0 0 0 0 0 1 0 0 0 44
45 1762444 0 0 0 0 0 0 0 0 1 0 0 45
46 1622044 0 0 0 0 0 0 0 0 0 1 0 46
47 1368955 0 0 0 0 0 0 0 0 0 0 1 47
48 1262973 0 0 0 0 0 0 0 0 0 0 0 48
49 1269530 1 0 0 0 0 0 0 0 0 0 0 49
50 1479279 0 1 0 0 0 0 0 0 0 0 0 50
51 1607819 0 0 1 0 0 0 0 0 0 0 0 51
52 1721466 0 0 0 1 0 0 0 0 0 0 0 52
53 1721766 0 0 0 0 1 0 0 0 0 0 0 53
54 1949843 0 0 0 0 0 1 0 0 0 0 0 54
55 1821326 0 0 0 0 0 0 1 0 0 0 0 55
56 1757802 0 0 0 0 0 0 0 1 0 0 0 56
57 1590367 0 0 0 0 0 0 0 0 1 0 0 57
58 1260647 0 0 0 0 0 0 0 0 0 1 0 58
59 1149235 0 0 0 0 0 0 0 0 0 0 1 59
60 1016367 0 0 0 0 0 0 0 0 0 0 0 60
61 1027885 1 0 0 0 0 0 0 0 0 0 0 61
62 1262159 0 1 0 0 0 0 0 0 0 0 0 62
63 1520854 0 0 1 0 0 0 0 0 0 0 0 63
64 1544144 0 0 0 1 0 0 0 0 0 0 0 64
65 1564709 0 0 0 0 1 0 0 0 0 0 0 65
66 1821776 0 0 0 0 0 1 0 0 0 0 0 66
67 1741365 0 0 0 0 0 0 1 0 0 0 0 67
68 1623386 0 0 0 0 0 0 0 1 0 0 0 68
69 1498658 0 0 0 0 0 0 0 0 1 0 0 69
70 1241822 0 0 0 0 0 0 0 0 0 1 0 70
71 1136029 0 0 0 0 0 0 0 0 0 0 1 71
72 1035030 0 0 0 0 0 0 0 0 0 0 0 72
73 1078521 1 0 0 0 0 0 0 0 0 0 0 73
74 1279431 0 1 0 0 0 0 0 0 0 0 0 74
75 1171023 0 0 1 0 0 0 0 0 0 0 0 75
76 1573377 0 0 0 1 0 0 0 0 0 0 0 76
77 1589514 0 0 0 0 1 0 0 0 0 0 0 77
78 1859878 0 0 0 0 0 1 0 0 0 0 0 78
79 1783191 0 0 0 0 0 0 1 0 0 0 0 79
80 1689849 0 0 0 0 0 0 0 1 0 0 0 80
81 1619868 0 0 0 0 0 0 0 0 1 0 0 81
82 1323443 0 0 0 0 0 0 0 0 0 1 0 82
83 1177481 0 0 0 0 0 0 0 0 0 0 1 83
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
990651 -18426 90722 247715 426767 482861
M6 M7 M8 M9 M10 M11
598496 704448 607099 521915 304467 101485
t
2068
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-222479 -78496 -5890 44604 294482
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 990650.8 52910.7 18.723 < 2e-16 ***
M1 -18426.0 65182.2 -0.283 0.778253
M2 90722.3 65162.1 1.392 0.168252
M3 247715.0 65146.3 3.802 0.000303 ***
M4 426766.5 65135.1 6.552 8.11e-09 ***
M5 482861.0 65128.4 7.414 2.18e-10 ***
M6 598495.5 65126.2 9.190 1.19e-13 ***
M7 704448.2 65128.4 10.816 < 2e-16 ***
M8 607098.6 65135.1 9.321 6.85e-14 ***
M9 521915.1 65146.3 8.011 1.74e-11 ***
M10 304467.3 65162.1 4.672 1.40e-05 ***
M11 101484.6 65182.2 1.557 0.123996
t 2068.5 540.7 3.826 0.000280 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 117100 on 70 degrees of freedom
Multiple R-squared: 0.8418, Adjusted R-squared: 0.8147
F-statistic: 31.04 on 12 and 70 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,] 4.753377e-03 9.506755e-03 0.9952466226
[2,] 5.916831e-04 1.183366e-03 0.9994083169
[3,] 2.305195e-04 4.610389e-04 0.9997694805
[4,] 3.201031e-05 6.402062e-05 0.9999679897
[5,] 1.406753e-05 2.813506e-05 0.9999859325
[6,] 1.968841e-06 3.937682e-06 0.9999980312
[7,] 4.621330e-07 9.242660e-07 0.9999995379
[8,] 7.047899e-08 1.409580e-07 0.9999999295
[9,] 4.926262e-08 9.852524e-08 0.9999999507
[10,] 7.493961e-09 1.498792e-08 0.9999999925
[11,] 2.632884e-09 5.265769e-09 0.9999999974
[12,] 8.003148e-10 1.600630e-09 0.9999999992
[13,] 1.196523e-09 2.393045e-09 0.9999999988
[14,] 4.920019e-10 9.840038e-10 0.9999999995
[15,] 2.784509e-10 5.569019e-10 0.9999999997
[16,] 6.808004e-11 1.361601e-10 0.9999999999
[17,] 1.254477e-10 2.508953e-10 0.9999999999
[18,] 2.504805e-11 5.009611e-11 1.0000000000
[19,] 4.837642e-12 9.675284e-12 1.0000000000
[20,] 9.810880e-13 1.962176e-12 1.0000000000
[21,] 3.124155e-13 6.248309e-13 1.0000000000
[22,] 1.153233e-13 2.306466e-13 1.0000000000
[23,] 4.897568e-13 9.795135e-13 1.0000000000
[24,] 1.228559e-12 2.457117e-12 1.0000000000
[25,] 9.611705e-13 1.922341e-12 1.0000000000
[26,] 4.976458e-13 9.952916e-13 1.0000000000
[27,] 1.368129e-10 2.736259e-10 0.9999999999
[28,] 1.746088e-10 3.492175e-10 0.9999999998
[29,] 1.522865e-10 3.045730e-10 0.9999999998
[30,] 3.679580e-10 7.359159e-10 0.9999999996
[31,] 1.203807e-08 2.407613e-08 0.9999999880
[32,] 6.866590e-08 1.373318e-07 0.9999999313
[33,] 6.132734e-07 1.226547e-06 0.9999993867
[34,] 1.159056e-05 2.318112e-05 0.9999884094
[35,] 1.505587e-02 3.011174e-02 0.9849441292
[36,] 1.136962e-01 2.273925e-01 0.8863037658
[37,] 1.601445e-01 3.202891e-01 0.8398554746
[38,] 1.719214e-01 3.438428e-01 0.8280785837
[39,] 3.247154e-01 6.494307e-01 0.6752846251
[40,] 3.711228e-01 7.422456e-01 0.6288771987
[41,] 4.247483e-01 8.494966e-01 0.5752516812
[42,] 5.084129e-01 9.831742e-01 0.4915871213
[43,] 7.156596e-01 5.686809e-01 0.2843404348
[44,] 7.005981e-01 5.988038e-01 0.2994018832
[45,] 6.774588e-01 6.450824e-01 0.3225411998
[46,] 6.265446e-01 7.469109e-01 0.3734554261
[47,] 5.186584e-01 9.626832e-01 0.4813416219
[48,] 9.997902e-01 4.195087e-04 0.0002097544
[49,] 9.993257e-01 1.348507e-03 0.0006742533
[50,] 9.983790e-01 3.242078e-03 0.0016210390
[51,] 9.946169e-01 1.076615e-02 0.0053830749
[52,] 9.843329e-01 3.133426e-02 0.0156671316
> postscript(file="/var/wessaorg/rcomp/tmp/1g9sr1323257416.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/2h3141323257416.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/3yf841323257416.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/450u61323257416.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/5wmpo1323257416.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 = 83
Frequency = 1
1 2 3 4 5 6
-52928.300 -97589.157 -111957.300 -93467.300 -65721.300 -190008.300
7 8 9 10 11 12
-13658.442 21875.700 8470.700 79510.986 12686.272 20603.346
13 14 15 16 17 18
-9879.104 -101951.962 -61630.104 -81674.104 -38878.104 -127658.104
19 20 21 22 23 24
27360.753 14094.896 43099.896 80554.181 24284.467 -2559.459
25 26 27 28 29 30
-8529.909 -95943.767 -36165.909 -5889.909 18848.091 -102056.909
31 32 33 34 35 36
26672.948 -1605.909 48614.091 101983.376 37676.662 11865.736
37 38 39 40 41 42
-9391.714 -96526.571 16098.286 -8554.714 33652.286 -34774.714
43 44 45 46 47 48
114805.143 109817.286 156796.286 231775.571 179600.857 173034.931
49 50 51 52 53 54
195949.481 294481.624 263960.481 196487.481 138624.481 248998.481
55 56 57 58 59 60
12460.338 44217.481 -40102.519 -154443.233 -64940.948 -98392.874
61 62 63 64 65 66
-70517.324 52539.819 152173.676 -5656.324 -43254.324 96109.676
67 68 69 70 71 72
-92322.467 -115020.324 -156633.324 -198090.038 -102968.753 -104551.679
73 74 75 76 77 78
-44703.129 44990.014 -222479.129 -1245.129 -43271.129 109389.871
79 80 81 82 83
-75318.272 -73379.129 -60245.129 -141290.843 -86338.558
> postscript(file="/var/wessaorg/rcomp/tmp/677fd1323257416.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 = 83
Frequency = 1
lag(myerror, k = 1) myerror
0 -52928.300 NA
1 -97589.157 -52928.300
2 -111957.300 -97589.157
3 -93467.300 -111957.300
4 -65721.300 -93467.300
5 -190008.300 -65721.300
6 -13658.442 -190008.300
7 21875.700 -13658.442
8 8470.700 21875.700
9 79510.986 8470.700
10 12686.272 79510.986
11 20603.346 12686.272
12 -9879.104 20603.346
13 -101951.962 -9879.104
14 -61630.104 -101951.962
15 -81674.104 -61630.104
16 -38878.104 -81674.104
17 -127658.104 -38878.104
18 27360.753 -127658.104
19 14094.896 27360.753
20 43099.896 14094.896
21 80554.181 43099.896
22 24284.467 80554.181
23 -2559.459 24284.467
24 -8529.909 -2559.459
25 -95943.767 -8529.909
26 -36165.909 -95943.767
27 -5889.909 -36165.909
28 18848.091 -5889.909
29 -102056.909 18848.091
30 26672.948 -102056.909
31 -1605.909 26672.948
32 48614.091 -1605.909
33 101983.376 48614.091
34 37676.662 101983.376
35 11865.736 37676.662
36 -9391.714 11865.736
37 -96526.571 -9391.714
38 16098.286 -96526.571
39 -8554.714 16098.286
40 33652.286 -8554.714
41 -34774.714 33652.286
42 114805.143 -34774.714
43 109817.286 114805.143
44 156796.286 109817.286
45 231775.571 156796.286
46 179600.857 231775.571
47 173034.931 179600.857
48 195949.481 173034.931
49 294481.624 195949.481
50 263960.481 294481.624
51 196487.481 263960.481
52 138624.481 196487.481
53 248998.481 138624.481
54 12460.338 248998.481
55 44217.481 12460.338
56 -40102.519 44217.481
57 -154443.233 -40102.519
58 -64940.948 -154443.233
59 -98392.874 -64940.948
60 -70517.324 -98392.874
61 52539.819 -70517.324
62 152173.676 52539.819
63 -5656.324 152173.676
64 -43254.324 -5656.324
65 96109.676 -43254.324
66 -92322.467 96109.676
67 -115020.324 -92322.467
68 -156633.324 -115020.324
69 -198090.038 -156633.324
70 -102968.753 -198090.038
71 -104551.679 -102968.753
72 -44703.129 -104551.679
73 44990.014 -44703.129
74 -222479.129 44990.014
75 -1245.129 -222479.129
76 -43271.129 -1245.129
77 109389.871 -43271.129
78 -75318.272 109389.871
79 -73379.129 -75318.272
80 -60245.129 -73379.129
81 -141290.843 -60245.129
82 -86338.558 -141290.843
83 NA -86338.558
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -97589.157 -52928.300
[2,] -111957.300 -97589.157
[3,] -93467.300 -111957.300
[4,] -65721.300 -93467.300
[5,] -190008.300 -65721.300
[6,] -13658.442 -190008.300
[7,] 21875.700 -13658.442
[8,] 8470.700 21875.700
[9,] 79510.986 8470.700
[10,] 12686.272 79510.986
[11,] 20603.346 12686.272
[12,] -9879.104 20603.346
[13,] -101951.962 -9879.104
[14,] -61630.104 -101951.962
[15,] -81674.104 -61630.104
[16,] -38878.104 -81674.104
[17,] -127658.104 -38878.104
[18,] 27360.753 -127658.104
[19,] 14094.896 27360.753
[20,] 43099.896 14094.896
[21,] 80554.181 43099.896
[22,] 24284.467 80554.181
[23,] -2559.459 24284.467
[24,] -8529.909 -2559.459
[25,] -95943.767 -8529.909
[26,] -36165.909 -95943.767
[27,] -5889.909 -36165.909
[28,] 18848.091 -5889.909
[29,] -102056.909 18848.091
[30,] 26672.948 -102056.909
[31,] -1605.909 26672.948
[32,] 48614.091 -1605.909
[33,] 101983.376 48614.091
[34,] 37676.662 101983.376
[35,] 11865.736 37676.662
[36,] -9391.714 11865.736
[37,] -96526.571 -9391.714
[38,] 16098.286 -96526.571
[39,] -8554.714 16098.286
[40,] 33652.286 -8554.714
[41,] -34774.714 33652.286
[42,] 114805.143 -34774.714
[43,] 109817.286 114805.143
[44,] 156796.286 109817.286
[45,] 231775.571 156796.286
[46,] 179600.857 231775.571
[47,] 173034.931 179600.857
[48,] 195949.481 173034.931
[49,] 294481.624 195949.481
[50,] 263960.481 294481.624
[51,] 196487.481 263960.481
[52,] 138624.481 196487.481
[53,] 248998.481 138624.481
[54,] 12460.338 248998.481
[55,] 44217.481 12460.338
[56,] -40102.519 44217.481
[57,] -154443.233 -40102.519
[58,] -64940.948 -154443.233
[59,] -98392.874 -64940.948
[60,] -70517.324 -98392.874
[61,] 52539.819 -70517.324
[62,] 152173.676 52539.819
[63,] -5656.324 152173.676
[64,] -43254.324 -5656.324
[65,] 96109.676 -43254.324
[66,] -92322.467 96109.676
[67,] -115020.324 -92322.467
[68,] -156633.324 -115020.324
[69,] -198090.038 -156633.324
[70,] -102968.753 -198090.038
[71,] -104551.679 -102968.753
[72,] -44703.129 -104551.679
[73,] 44990.014 -44703.129
[74,] -222479.129 44990.014
[75,] -1245.129 -222479.129
[76,] -43271.129 -1245.129
[77,] 109389.871 -43271.129
[78,] -75318.272 109389.871
[79,] -73379.129 -75318.272
[80,] -60245.129 -73379.129
[81,] -141290.843 -60245.129
[82,] -86338.558 -141290.843
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -97589.157 -52928.300
2 -111957.300 -97589.157
3 -93467.300 -111957.300
4 -65721.300 -93467.300
5 -190008.300 -65721.300
6 -13658.442 -190008.300
7 21875.700 -13658.442
8 8470.700 21875.700
9 79510.986 8470.700
10 12686.272 79510.986
11 20603.346 12686.272
12 -9879.104 20603.346
13 -101951.962 -9879.104
14 -61630.104 -101951.962
15 -81674.104 -61630.104
16 -38878.104 -81674.104
17 -127658.104 -38878.104
18 27360.753 -127658.104
19 14094.896 27360.753
20 43099.896 14094.896
21 80554.181 43099.896
22 24284.467 80554.181
23 -2559.459 24284.467
24 -8529.909 -2559.459
25 -95943.767 -8529.909
26 -36165.909 -95943.767
27 -5889.909 -36165.909
28 18848.091 -5889.909
29 -102056.909 18848.091
30 26672.948 -102056.909
31 -1605.909 26672.948
32 48614.091 -1605.909
33 101983.376 48614.091
34 37676.662 101983.376
35 11865.736 37676.662
36 -9391.714 11865.736
37 -96526.571 -9391.714
38 16098.286 -96526.571
39 -8554.714 16098.286
40 33652.286 -8554.714
41 -34774.714 33652.286
42 114805.143 -34774.714
43 109817.286 114805.143
44 156796.286 109817.286
45 231775.571 156796.286
46 179600.857 231775.571
47 173034.931 179600.857
48 195949.481 173034.931
49 294481.624 195949.481
50 263960.481 294481.624
51 196487.481 263960.481
52 138624.481 196487.481
53 248998.481 138624.481
54 12460.338 248998.481
55 44217.481 12460.338
56 -40102.519 44217.481
57 -154443.233 -40102.519
58 -64940.948 -154443.233
59 -98392.874 -64940.948
60 -70517.324 -98392.874
61 52539.819 -70517.324
62 152173.676 52539.819
63 -5656.324 152173.676
64 -43254.324 -5656.324
65 96109.676 -43254.324
66 -92322.467 96109.676
67 -115020.324 -92322.467
68 -156633.324 -115020.324
69 -198090.038 -156633.324
70 -102968.753 -198090.038
71 -104551.679 -102968.753
72 -44703.129 -104551.679
73 44990.014 -44703.129
74 -222479.129 44990.014
75 -1245.129 -222479.129
76 -43271.129 -1245.129
77 109389.871 -43271.129
78 -75318.272 109389.871
79 -73379.129 -75318.272
80 -60245.129 -73379.129
81 -141290.843 -60245.129
82 -86338.558 -141290.843
> 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/7th2v1323257416.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/8pwg61323257416.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/969yo1323257416.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/10d0sb1323257416.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/11f1gl1323257416.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/122bry1323257416.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/13xgm41323257416.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/14n2ih1323257416.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/15kkd91323257416.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/16gtgn1323257416.tab")
+ }
>
> try(system("convert tmp/1g9sr1323257416.ps tmp/1g9sr1323257416.png",intern=TRUE))
character(0)
> try(system("convert tmp/2h3141323257416.ps tmp/2h3141323257416.png",intern=TRUE))
character(0)
> try(system("convert tmp/3yf841323257416.ps tmp/3yf841323257416.png",intern=TRUE))
character(0)
> try(system("convert tmp/450u61323257416.ps tmp/450u61323257416.png",intern=TRUE))
character(0)
> try(system("convert tmp/5wmpo1323257416.ps tmp/5wmpo1323257416.png",intern=TRUE))
character(0)
> try(system("convert tmp/677fd1323257416.ps tmp/677fd1323257416.png",intern=TRUE))
character(0)
> try(system("convert tmp/7th2v1323257416.ps tmp/7th2v1323257416.png",intern=TRUE))
character(0)
> try(system("convert tmp/8pwg61323257416.ps tmp/8pwg61323257416.png",intern=TRUE))
character(0)
> try(system("convert tmp/969yo1323257416.ps tmp/969yo1323257416.png",intern=TRUE))
character(0)
> try(system("convert tmp/10d0sb1323257416.ps tmp/10d0sb1323257416.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.336 0.463 3.831