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(6827,6178,7084,8162,8462,9644,10466,10748,9963,8194,6848,7027,7269,6775,7819,8371,9069,10248,11030,10882,10333,9109,7685,7602,8350,7829,8829,9948,10638,11253,11424,11391,10665,9396,7775,7933,8186,7444,8484,9864,10252,12282,11637,11577,12417,9637,8094,9280,8334,7899,9994,10078,10801,12950,12222,12246,13281,10366,8730,9614,8639,8772,10894,10455,11179,10588,10794,12770,13812,10857,9290,10925,9491,8919,11607,8852,12537,14759,13667,13731,15110,12185,10645,12161,10840,10436,13589,13402,13103,14933,14147,14057,16234,12389,11595,12772),dim=c(1,96),dimnames=list(c('#Miles'),1:96))
> y <- array(NA,dim=c(1,96),dimnames=list(c('#Miles'),1:96))
> 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 = 'Do not include Seasonal 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
#Miles t
1 6827 1
2 6178 2
3 7084 3
4 8162 4
5 8462 5
6 9644 6
7 10466 7
8 10748 8
9 9963 9
10 8194 10
11 6848 11
12 7027 12
13 7269 13
14 6775 14
15 7819 15
16 8371 16
17 9069 17
18 10248 18
19 11030 19
20 10882 20
21 10333 21
22 9109 22
23 7685 23
24 7602 24
25 8350 25
26 7829 26
27 8829 27
28 9948 28
29 10638 29
30 11253 30
31 11424 31
32 11391 32
33 10665 33
34 9396 34
35 7775 35
36 7933 36
37 8186 37
38 7444 38
39 8484 39
40 9864 40
41 10252 41
42 12282 42
43 11637 43
44 11577 44
45 12417 45
46 9637 46
47 8094 47
48 9280 48
49 8334 49
50 7899 50
51 9994 51
52 10078 52
53 10801 53
54 12950 54
55 12222 55
56 12246 56
57 13281 57
58 10366 58
59 8730 59
60 9614 60
61 8639 61
62 8772 62
63 10894 63
64 10455 64
65 11179 65
66 10588 66
67 10794 67
68 12770 68
69 13812 69
70 10857 70
71 9290 71
72 10925 72
73 9491 73
74 8919 74
75 11607 75
76 8852 76
77 12537 77
78 14759 78
79 13667 79
80 13731 80
81 15110 81
82 12185 82
83 10645 83
84 12161 84
85 10840 85
86 10436 86
87 13589 87
88 13402 88
89 13103 89
90 14933 90
91 14147 91
92 14057 92
93 16234 93
94 12389 94
95 11595 95
96 12772 96
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) t
7747.2 54.4
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3029.8 -1342.8 -199.4 1449.9 3427.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7747.174 330.542 23.438 < 2e-16 ***
t 54.403 5.917 9.194 9.53e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1607 on 94 degrees of freedom
Multiple R-squared: 0.4735, Adjusted R-squared: 0.4679
F-statistic: 84.52 on 1 and 94 DF, p-value: 9.528e-15
> 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.043926743 0.087853485 0.9560733
[2,] 0.018869820 0.037739640 0.9811302
[3,] 0.007265390 0.014530781 0.9927346
[4,] 0.002143325 0.004286649 0.9978567
[5,] 0.006564074 0.013128148 0.9934359
[6,] 0.126676651 0.253353302 0.8733233
[7,] 0.432658608 0.865317216 0.5673414
[8,] 0.516758451 0.966483099 0.4832415
[9,] 0.498472142 0.996944285 0.5015279
[10,] 0.494671739 0.989343478 0.5053283
[11,] 0.411481923 0.822963847 0.5885181
[12,] 0.331935249 0.663870498 0.6680648
[13,] 0.276263899 0.552527798 0.7237361
[14,] 0.285473298 0.570946596 0.7145267
[15,] 0.334896058 0.669792116 0.6651039
[16,] 0.334394436 0.668788872 0.6656056
[17,] 0.288071886 0.576143773 0.7119281
[18,] 0.240266479 0.480532959 0.7597335
[19,] 0.268559154 0.537118307 0.7314408
[20,] 0.283219941 0.566439882 0.7167801
[21,] 0.243460829 0.486921658 0.7565392
[22,] 0.226787673 0.453575346 0.7732123
[23,] 0.179394041 0.358788082 0.8206060
[24,] 0.147779202 0.295558404 0.8522208
[25,] 0.137354013 0.274708025 0.8626460
[26,] 0.148371377 0.296742755 0.8516286
[27,] 0.160660966 0.321321933 0.8393390
[28,] 0.165579801 0.331159602 0.8344202
[29,] 0.142451150 0.284902300 0.8575489
[30,] 0.119943762 0.239887523 0.8800562
[31,] 0.156093900 0.312187799 0.8439061
[32,] 0.174385075 0.348770149 0.8256149
[33,] 0.172821517 0.345643034 0.8271785
[34,] 0.211057368 0.422114737 0.7889426
[35,] 0.188885268 0.377770535 0.8111147
[36,] 0.151143509 0.302287019 0.8488565
[37,] 0.122021034 0.244042069 0.8779790
[38,] 0.172192766 0.344385533 0.8278072
[39,] 0.180595416 0.361190832 0.8194046
[40,] 0.183720168 0.367440336 0.8162798
[41,] 0.244545581 0.489091161 0.7554544
[42,] 0.209975642 0.419951285 0.7900244
[43,] 0.237707431 0.475414861 0.7622926
[44,] 0.206481667 0.412963334 0.7935183
[45,] 0.214739341 0.429478683 0.7852607
[46,] 0.253063131 0.506126261 0.7469369
[47,] 0.208210848 0.416421697 0.7917892
[48,] 0.168305613 0.336611227 0.8316944
[49,] 0.137353067 0.274706134 0.8626469
[50,] 0.201066291 0.402132581 0.7989337
[51,] 0.220031396 0.440062792 0.7799686
[52,] 0.244367717 0.488735434 0.7556323
[53,] 0.392582948 0.785165896 0.6074171
[54,] 0.348069552 0.696139103 0.6519304
[55,] 0.352173847 0.704347694 0.6478262
[56,] 0.312976430 0.625952860 0.6870236
[57,] 0.326918586 0.653837171 0.6730814
[58,] 0.337675373 0.675350746 0.6623246
[59,] 0.284424958 0.568849916 0.7155750
[60,] 0.234911783 0.469823566 0.7650882
[61,] 0.192263313 0.384526626 0.8077367
[62,] 0.152742211 0.305484422 0.8472578
[63,] 0.118305510 0.236611019 0.8816945
[64,] 0.126155960 0.252311921 0.8738440
[65,] 0.222503771 0.445007543 0.7774962
[66,] 0.178601851 0.357203701 0.8213981
[67,] 0.177427448 0.354854896 0.8225726
[68,] 0.137299265 0.274598530 0.8627007
[69,] 0.141342132 0.282684264 0.8586579
[70,] 0.214484455 0.428968911 0.7855155
[71,] 0.170497204 0.340994408 0.8295028
[72,] 0.381553803 0.763107605 0.6184462
[73,] 0.328040745 0.656081490 0.6719593
[74,] 0.395091421 0.790182843 0.6049086
[75,] 0.367896848 0.735793697 0.6321032
[76,] 0.352117139 0.704234278 0.6478829
[77,] 0.586832144 0.826335713 0.4131679
[78,] 0.508968325 0.982063351 0.4910317
[79,] 0.471277686 0.942555371 0.5287223
[80,] 0.378718282 0.757436563 0.6212817
[81,] 0.390200780 0.780401560 0.6097992
[82,] 0.679282615 0.641434771 0.3207174
[83,] 0.599489370 0.801021261 0.4005106
[84,] 0.543404848 0.913190304 0.4565952
[85,] 0.611258087 0.777483827 0.3887419
[86,] 0.483127246 0.966254491 0.5168728
[87,] 0.379648283 0.759296566 0.6203517
> postscript(file="/var/wessaorg/rcomp/tmp/1pch51322313564.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/2ddiw1322313564.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/3exz81322313564.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/4u2uz1322313564.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/5g4un1322313564.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 = 96
Frequency = 1
1 2 3 4 5 6
-974.57753 -1677.98095 -826.38436 197.21223 442.80882 1570.40541
7 8 9 10 11 12
2338.00200 2565.59858 1726.19517 -97.20824 -1497.61165 -1373.01506
13 14 15 16 17 18
-1185.41847 -1733.82188 -744.22530 -246.62871 396.96788 1521.56447
19 20 21 22 23 24
2249.16106 2046.75765 1443.35423 164.95082 -1313.45259 -1450.85600
25 26 27 28 29 30
-757.25941 -1332.66282 -387.06623 677.53035 1313.12694 1873.72353
31 32 33 34 35 36
1990.32012 1902.91671 1122.51330 -200.89012 -1876.29353 -1772.69694
37 38 39 40 41 42
-1574.10035 -2370.50376 -1384.90717 -59.31059 274.28600 2249.88259
43 44 45 46 47 48
1550.47918 1436.07577 2221.67236 -612.73105 -2210.13447 -1078.53788
49 50 51 52 53 54
-2078.94129 -2568.34470 -527.74811 -498.15152 170.44506 2265.04165
55 56 57 58 59 60
1482.63824 1452.23483 2432.83142 -536.57199 -2226.97540 -1397.37882
61 62 63 64 65 66
-2426.78223 -2348.18564 -280.58905 -773.99246 -104.39587 -749.79929
67 68 69 70 71 72
-598.20270 1323.39389 2310.99048 -698.41293 -2319.81634 -739.21975
73 74 75 76 77 78
-2227.62317 -2854.02658 -220.42999 -3029.83340 600.76319 2768.35978
79 80 81 82 83 84
1621.95636 1631.55295 2956.14954 -23.25387 -1617.65728 -156.06069
85 86 87 88 89 90
-1531.46411 -1989.86752 1108.72907 867.32566 513.92225 2289.51884
91 92 93 94 95 96
1449.11543 1304.71201 3427.30860 -472.09481 -1320.49822 -197.90163
> postscript(file="/var/wessaorg/rcomp/tmp/617i81322313564.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 = 96
Frequency = 1
lag(myerror, k = 1) myerror
0 -974.57753 NA
1 -1677.98095 -974.57753
2 -826.38436 -1677.98095
3 197.21223 -826.38436
4 442.80882 197.21223
5 1570.40541 442.80882
6 2338.00200 1570.40541
7 2565.59858 2338.00200
8 1726.19517 2565.59858
9 -97.20824 1726.19517
10 -1497.61165 -97.20824
11 -1373.01506 -1497.61165
12 -1185.41847 -1373.01506
13 -1733.82188 -1185.41847
14 -744.22530 -1733.82188
15 -246.62871 -744.22530
16 396.96788 -246.62871
17 1521.56447 396.96788
18 2249.16106 1521.56447
19 2046.75765 2249.16106
20 1443.35423 2046.75765
21 164.95082 1443.35423
22 -1313.45259 164.95082
23 -1450.85600 -1313.45259
24 -757.25941 -1450.85600
25 -1332.66282 -757.25941
26 -387.06623 -1332.66282
27 677.53035 -387.06623
28 1313.12694 677.53035
29 1873.72353 1313.12694
30 1990.32012 1873.72353
31 1902.91671 1990.32012
32 1122.51330 1902.91671
33 -200.89012 1122.51330
34 -1876.29353 -200.89012
35 -1772.69694 -1876.29353
36 -1574.10035 -1772.69694
37 -2370.50376 -1574.10035
38 -1384.90717 -2370.50376
39 -59.31059 -1384.90717
40 274.28600 -59.31059
41 2249.88259 274.28600
42 1550.47918 2249.88259
43 1436.07577 1550.47918
44 2221.67236 1436.07577
45 -612.73105 2221.67236
46 -2210.13447 -612.73105
47 -1078.53788 -2210.13447
48 -2078.94129 -1078.53788
49 -2568.34470 -2078.94129
50 -527.74811 -2568.34470
51 -498.15152 -527.74811
52 170.44506 -498.15152
53 2265.04165 170.44506
54 1482.63824 2265.04165
55 1452.23483 1482.63824
56 2432.83142 1452.23483
57 -536.57199 2432.83142
58 -2226.97540 -536.57199
59 -1397.37882 -2226.97540
60 -2426.78223 -1397.37882
61 -2348.18564 -2426.78223
62 -280.58905 -2348.18564
63 -773.99246 -280.58905
64 -104.39587 -773.99246
65 -749.79929 -104.39587
66 -598.20270 -749.79929
67 1323.39389 -598.20270
68 2310.99048 1323.39389
69 -698.41293 2310.99048
70 -2319.81634 -698.41293
71 -739.21975 -2319.81634
72 -2227.62317 -739.21975
73 -2854.02658 -2227.62317
74 -220.42999 -2854.02658
75 -3029.83340 -220.42999
76 600.76319 -3029.83340
77 2768.35978 600.76319
78 1621.95636 2768.35978
79 1631.55295 1621.95636
80 2956.14954 1631.55295
81 -23.25387 2956.14954
82 -1617.65728 -23.25387
83 -156.06069 -1617.65728
84 -1531.46411 -156.06069
85 -1989.86752 -1531.46411
86 1108.72907 -1989.86752
87 867.32566 1108.72907
88 513.92225 867.32566
89 2289.51884 513.92225
90 1449.11543 2289.51884
91 1304.71201 1449.11543
92 3427.30860 1304.71201
93 -472.09481 3427.30860
94 -1320.49822 -472.09481
95 -197.90163 -1320.49822
96 NA -197.90163
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1677.98095 -974.57753
[2,] -826.38436 -1677.98095
[3,] 197.21223 -826.38436
[4,] 442.80882 197.21223
[5,] 1570.40541 442.80882
[6,] 2338.00200 1570.40541
[7,] 2565.59858 2338.00200
[8,] 1726.19517 2565.59858
[9,] -97.20824 1726.19517
[10,] -1497.61165 -97.20824
[11,] -1373.01506 -1497.61165
[12,] -1185.41847 -1373.01506
[13,] -1733.82188 -1185.41847
[14,] -744.22530 -1733.82188
[15,] -246.62871 -744.22530
[16,] 396.96788 -246.62871
[17,] 1521.56447 396.96788
[18,] 2249.16106 1521.56447
[19,] 2046.75765 2249.16106
[20,] 1443.35423 2046.75765
[21,] 164.95082 1443.35423
[22,] -1313.45259 164.95082
[23,] -1450.85600 -1313.45259
[24,] -757.25941 -1450.85600
[25,] -1332.66282 -757.25941
[26,] -387.06623 -1332.66282
[27,] 677.53035 -387.06623
[28,] 1313.12694 677.53035
[29,] 1873.72353 1313.12694
[30,] 1990.32012 1873.72353
[31,] 1902.91671 1990.32012
[32,] 1122.51330 1902.91671
[33,] -200.89012 1122.51330
[34,] -1876.29353 -200.89012
[35,] -1772.69694 -1876.29353
[36,] -1574.10035 -1772.69694
[37,] -2370.50376 -1574.10035
[38,] -1384.90717 -2370.50376
[39,] -59.31059 -1384.90717
[40,] 274.28600 -59.31059
[41,] 2249.88259 274.28600
[42,] 1550.47918 2249.88259
[43,] 1436.07577 1550.47918
[44,] 2221.67236 1436.07577
[45,] -612.73105 2221.67236
[46,] -2210.13447 -612.73105
[47,] -1078.53788 -2210.13447
[48,] -2078.94129 -1078.53788
[49,] -2568.34470 -2078.94129
[50,] -527.74811 -2568.34470
[51,] -498.15152 -527.74811
[52,] 170.44506 -498.15152
[53,] 2265.04165 170.44506
[54,] 1482.63824 2265.04165
[55,] 1452.23483 1482.63824
[56,] 2432.83142 1452.23483
[57,] -536.57199 2432.83142
[58,] -2226.97540 -536.57199
[59,] -1397.37882 -2226.97540
[60,] -2426.78223 -1397.37882
[61,] -2348.18564 -2426.78223
[62,] -280.58905 -2348.18564
[63,] -773.99246 -280.58905
[64,] -104.39587 -773.99246
[65,] -749.79929 -104.39587
[66,] -598.20270 -749.79929
[67,] 1323.39389 -598.20270
[68,] 2310.99048 1323.39389
[69,] -698.41293 2310.99048
[70,] -2319.81634 -698.41293
[71,] -739.21975 -2319.81634
[72,] -2227.62317 -739.21975
[73,] -2854.02658 -2227.62317
[74,] -220.42999 -2854.02658
[75,] -3029.83340 -220.42999
[76,] 600.76319 -3029.83340
[77,] 2768.35978 600.76319
[78,] 1621.95636 2768.35978
[79,] 1631.55295 1621.95636
[80,] 2956.14954 1631.55295
[81,] -23.25387 2956.14954
[82,] -1617.65728 -23.25387
[83,] -156.06069 -1617.65728
[84,] -1531.46411 -156.06069
[85,] -1989.86752 -1531.46411
[86,] 1108.72907 -1989.86752
[87,] 867.32566 1108.72907
[88,] 513.92225 867.32566
[89,] 2289.51884 513.92225
[90,] 1449.11543 2289.51884
[91,] 1304.71201 1449.11543
[92,] 3427.30860 1304.71201
[93,] -472.09481 3427.30860
[94,] -1320.49822 -472.09481
[95,] -197.90163 -1320.49822
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1677.98095 -974.57753
2 -826.38436 -1677.98095
3 197.21223 -826.38436
4 442.80882 197.21223
5 1570.40541 442.80882
6 2338.00200 1570.40541
7 2565.59858 2338.00200
8 1726.19517 2565.59858
9 -97.20824 1726.19517
10 -1497.61165 -97.20824
11 -1373.01506 -1497.61165
12 -1185.41847 -1373.01506
13 -1733.82188 -1185.41847
14 -744.22530 -1733.82188
15 -246.62871 -744.22530
16 396.96788 -246.62871
17 1521.56447 396.96788
18 2249.16106 1521.56447
19 2046.75765 2249.16106
20 1443.35423 2046.75765
21 164.95082 1443.35423
22 -1313.45259 164.95082
23 -1450.85600 -1313.45259
24 -757.25941 -1450.85600
25 -1332.66282 -757.25941
26 -387.06623 -1332.66282
27 677.53035 -387.06623
28 1313.12694 677.53035
29 1873.72353 1313.12694
30 1990.32012 1873.72353
31 1902.91671 1990.32012
32 1122.51330 1902.91671
33 -200.89012 1122.51330
34 -1876.29353 -200.89012
35 -1772.69694 -1876.29353
36 -1574.10035 -1772.69694
37 -2370.50376 -1574.10035
38 -1384.90717 -2370.50376
39 -59.31059 -1384.90717
40 274.28600 -59.31059
41 2249.88259 274.28600
42 1550.47918 2249.88259
43 1436.07577 1550.47918
44 2221.67236 1436.07577
45 -612.73105 2221.67236
46 -2210.13447 -612.73105
47 -1078.53788 -2210.13447
48 -2078.94129 -1078.53788
49 -2568.34470 -2078.94129
50 -527.74811 -2568.34470
51 -498.15152 -527.74811
52 170.44506 -498.15152
53 2265.04165 170.44506
54 1482.63824 2265.04165
55 1452.23483 1482.63824
56 2432.83142 1452.23483
57 -536.57199 2432.83142
58 -2226.97540 -536.57199
59 -1397.37882 -2226.97540
60 -2426.78223 -1397.37882
61 -2348.18564 -2426.78223
62 -280.58905 -2348.18564
63 -773.99246 -280.58905
64 -104.39587 -773.99246
65 -749.79929 -104.39587
66 -598.20270 -749.79929
67 1323.39389 -598.20270
68 2310.99048 1323.39389
69 -698.41293 2310.99048
70 -2319.81634 -698.41293
71 -739.21975 -2319.81634
72 -2227.62317 -739.21975
73 -2854.02658 -2227.62317
74 -220.42999 -2854.02658
75 -3029.83340 -220.42999
76 600.76319 -3029.83340
77 2768.35978 600.76319
78 1621.95636 2768.35978
79 1631.55295 1621.95636
80 2956.14954 1631.55295
81 -23.25387 2956.14954
82 -1617.65728 -23.25387
83 -156.06069 -1617.65728
84 -1531.46411 -156.06069
85 -1989.86752 -1531.46411
86 1108.72907 -1989.86752
87 867.32566 1108.72907
88 513.92225 867.32566
89 2289.51884 513.92225
90 1449.11543 2289.51884
91 1304.71201 1449.11543
92 3427.30860 1304.71201
93 -472.09481 3427.30860
94 -1320.49822 -472.09481
95 -197.90163 -1320.49822
> 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/7rsk61322313564.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/8klga1322313564.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/9yyfb1322313564.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/10cx8q1322313564.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/11v8yp1322313564.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/12bspr1322313564.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/13xkfz1322313564.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/14rx8j1322313564.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/15qspa1322313564.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/16e6kd1322313564.tab")
+ }
>
> try(system("convert tmp/1pch51322313564.ps tmp/1pch51322313564.png",intern=TRUE))
character(0)
> try(system("convert tmp/2ddiw1322313564.ps tmp/2ddiw1322313564.png",intern=TRUE))
character(0)
> try(system("convert tmp/3exz81322313564.ps tmp/3exz81322313564.png",intern=TRUE))
character(0)
> try(system("convert tmp/4u2uz1322313564.ps tmp/4u2uz1322313564.png",intern=TRUE))
character(0)
> try(system("convert tmp/5g4un1322313564.ps tmp/5g4un1322313564.png",intern=TRUE))
character(0)
> try(system("convert tmp/617i81322313564.ps tmp/617i81322313564.png",intern=TRUE))
character(0)
> try(system("convert tmp/7rsk61322313564.ps tmp/7rsk61322313564.png",intern=TRUE))
character(0)
> try(system("convert tmp/8klga1322313564.ps tmp/8klga1322313564.png",intern=TRUE))
character(0)
> try(system("convert tmp/9yyfb1322313564.ps tmp/9yyfb1322313564.png",intern=TRUE))
character(0)
> try(system("convert tmp/10cx8q1322313564.ps tmp/10cx8q1322313564.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.755 0.493 4.329