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(149657,0,142773,0,133639,0,128332,0,120297,0,118632,0,155276,0,169316,0,167395,0,157939,0,149601,0,146310,0,141579,0,136473,0,129818,0,124226,0,116428,0,116440,0,147747,0,160069,0,163129,0,151108,0,141481,0,139174,0,134066,0,130104,0,123090,0,116598,0,109627,0,105428,0,137272,0,159836,0,155283,0,141514,0,131852,0,130691,0,128461,0,123066,0,117599,0,111599,0,105395,0,102334,0,131305,0,149033,0,144954,0,132404,0,122104,0,118755,0,116222,1,110924,1,103753,1,99983,1,93302,1,91496,1,119321,1,139261,1,133739,1,123913,1,113438,1,109416,1,109406,1,105645,1,101328,1,97686,1,93093,1,91382,1,122257,1,139183,1,139887,1,131822,1,116805,1,113706,1,113012,1,110452,1,107005,1,102841,1,98173,1,98181,1,137277,1,147579,1,146571,1,138920,1,130340,1,128140,1,127059,1,122860,1,117702,1,113537,1,108366,1,111078,1,150739,1,159129,1,157928,1,147768,1,137507,1,136919,1),dim=c(2,96),dimnames=list(c('Y','X'),1:96))
> y <- array(NA,dim=c(2,96),dimnames=list(c('Y','X'),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 = '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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 149657 0 1 0 0 0 0 0 0 0 0 0 0 1
2 142773 0 0 1 0 0 0 0 0 0 0 0 0 2
3 133639 0 0 0 1 0 0 0 0 0 0 0 0 3
4 128332 0 0 0 0 1 0 0 0 0 0 0 0 4
5 120297 0 0 0 0 0 1 0 0 0 0 0 0 5
6 118632 0 0 0 0 0 0 1 0 0 0 0 0 6
7 155276 0 0 0 0 0 0 0 1 0 0 0 0 7
8 169316 0 0 0 0 0 0 0 0 1 0 0 0 8
9 167395 0 0 0 0 0 0 0 0 0 1 0 0 9
10 157939 0 0 0 0 0 0 0 0 0 0 1 0 10
11 149601 0 0 0 0 0 0 0 0 0 0 0 1 11
12 146310 0 0 0 0 0 0 0 0 0 0 0 0 12
13 141579 0 1 0 0 0 0 0 0 0 0 0 0 13
14 136473 0 0 1 0 0 0 0 0 0 0 0 0 14
15 129818 0 0 0 1 0 0 0 0 0 0 0 0 15
16 124226 0 0 0 0 1 0 0 0 0 0 0 0 16
17 116428 0 0 0 0 0 1 0 0 0 0 0 0 17
18 116440 0 0 0 0 0 0 1 0 0 0 0 0 18
19 147747 0 0 0 0 0 0 0 1 0 0 0 0 19
20 160069 0 0 0 0 0 0 0 0 1 0 0 0 20
21 163129 0 0 0 0 0 0 0 0 0 1 0 0 21
22 151108 0 0 0 0 0 0 0 0 0 0 1 0 22
23 141481 0 0 0 0 0 0 0 0 0 0 0 1 23
24 139174 0 0 0 0 0 0 0 0 0 0 0 0 24
25 134066 0 1 0 0 0 0 0 0 0 0 0 0 25
26 130104 0 0 1 0 0 0 0 0 0 0 0 0 26
27 123090 0 0 0 1 0 0 0 0 0 0 0 0 27
28 116598 0 0 0 0 1 0 0 0 0 0 0 0 28
29 109627 0 0 0 0 0 1 0 0 0 0 0 0 29
30 105428 0 0 0 0 0 0 1 0 0 0 0 0 30
31 137272 0 0 0 0 0 0 0 1 0 0 0 0 31
32 159836 0 0 0 0 0 0 0 0 1 0 0 0 32
33 155283 0 0 0 0 0 0 0 0 0 1 0 0 33
34 141514 0 0 0 0 0 0 0 0 0 0 1 0 34
35 131852 0 0 0 0 0 0 0 0 0 0 0 1 35
36 130691 0 0 0 0 0 0 0 0 0 0 0 0 36
37 128461 0 1 0 0 0 0 0 0 0 0 0 0 37
38 123066 0 0 1 0 0 0 0 0 0 0 0 0 38
39 117599 0 0 0 1 0 0 0 0 0 0 0 0 39
40 111599 0 0 0 0 1 0 0 0 0 0 0 0 40
41 105395 0 0 0 0 0 1 0 0 0 0 0 0 41
42 102334 0 0 0 0 0 0 1 0 0 0 0 0 42
43 131305 0 0 0 0 0 0 0 1 0 0 0 0 43
44 149033 0 0 0 0 0 0 0 0 1 0 0 0 44
45 144954 0 0 0 0 0 0 0 0 0 1 0 0 45
46 132404 0 0 0 0 0 0 0 0 0 0 1 0 46
47 122104 0 0 0 0 0 0 0 0 0 0 0 1 47
48 118755 0 0 0 0 0 0 0 0 0 0 0 0 48
49 116222 1 1 0 0 0 0 0 0 0 0 0 0 49
50 110924 1 0 1 0 0 0 0 0 0 0 0 0 50
51 103753 1 0 0 1 0 0 0 0 0 0 0 0 51
52 99983 1 0 0 0 1 0 0 0 0 0 0 0 52
53 93302 1 0 0 0 0 1 0 0 0 0 0 0 53
54 91496 1 0 0 0 0 0 1 0 0 0 0 0 54
55 119321 1 0 0 0 0 0 0 1 0 0 0 0 55
56 139261 1 0 0 0 0 0 0 0 1 0 0 0 56
57 133739 1 0 0 0 0 0 0 0 0 1 0 0 57
58 123913 1 0 0 0 0 0 0 0 0 0 1 0 58
59 113438 1 0 0 0 0 0 0 0 0 0 0 1 59
60 109416 1 0 0 0 0 0 0 0 0 0 0 0 60
61 109406 1 1 0 0 0 0 0 0 0 0 0 0 61
62 105645 1 0 1 0 0 0 0 0 0 0 0 0 62
63 101328 1 0 0 1 0 0 0 0 0 0 0 0 63
64 97686 1 0 0 0 1 0 0 0 0 0 0 0 64
65 93093 1 0 0 0 0 1 0 0 0 0 0 0 65
66 91382 1 0 0 0 0 0 1 0 0 0 0 0 66
67 122257 1 0 0 0 0 0 0 1 0 0 0 0 67
68 139183 1 0 0 0 0 0 0 0 1 0 0 0 68
69 139887 1 0 0 0 0 0 0 0 0 1 0 0 69
70 131822 1 0 0 0 0 0 0 0 0 0 1 0 70
71 116805 1 0 0 0 0 0 0 0 0 0 0 1 71
72 113706 1 0 0 0 0 0 0 0 0 0 0 0 72
73 113012 1 1 0 0 0 0 0 0 0 0 0 0 73
74 110452 1 0 1 0 0 0 0 0 0 0 0 0 74
75 107005 1 0 0 1 0 0 0 0 0 0 0 0 75
76 102841 1 0 0 0 1 0 0 0 0 0 0 0 76
77 98173 1 0 0 0 0 1 0 0 0 0 0 0 77
78 98181 1 0 0 0 0 0 1 0 0 0 0 0 78
79 137277 1 0 0 0 0 0 0 1 0 0 0 0 79
80 147579 1 0 0 0 0 0 0 0 1 0 0 0 80
81 146571 1 0 0 0 0 0 0 0 0 1 0 0 81
82 138920 1 0 0 0 0 0 0 0 0 0 1 0 82
83 130340 1 0 0 0 0 0 0 0 0 0 0 1 83
84 128140 1 0 0 0 0 0 0 0 0 0 0 0 84
85 127059 1 1 0 0 0 0 0 0 0 0 0 0 85
86 122860 1 0 1 0 0 0 0 0 0 0 0 0 86
87 117702 1 0 0 1 0 0 0 0 0 0 0 0 87
88 113537 1 0 0 0 1 0 0 0 0 0 0 0 88
89 108366 1 0 0 0 0 1 0 0 0 0 0 0 89
90 111078 1 0 0 0 0 0 1 0 0 0 0 0 90
91 150739 1 0 0 0 0 0 0 1 0 0 0 0 91
92 159129 1 0 0 0 0 0 0 0 1 0 0 0 92
93 157928 1 0 0 0 0 0 0 0 0 1 0 0 93
94 147768 1 0 0 0 0 0 0 0 0 0 1 0 94
95 137507 1 0 0 0 0 0 0 0 0 0 0 1 95
96 136919 1 0 0 0 0 0 0 0 0 0 0 0 96
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
135682.42 -13898.76 -628.08 -5258.08 -11287.82 -16163.69
M5 M6 M7 M8 M9 M10
-22413.18 -23611.30 9682.21 24974.34 23174.98 12753.36
M11 t
2486.49 -15.63
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-16177 -6539 -1237 7137 20696
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 135682.42 3993.16 33.979 < 2e-16 ***
X -13898.76 3857.76 -3.603 0.000538 ***
M1 -628.08 4675.03 -0.134 0.893456
M2 -5258.08 4663.97 -1.127 0.262869
M3 -11287.82 4653.93 -2.425 0.017486 *
M4 -16163.69 4644.93 -3.480 0.000806 ***
M5 -22413.18 4636.98 -4.834 6.17e-06 ***
M6 -23611.30 4630.08 -5.100 2.15e-06 ***
M7 9682.21 4624.23 2.094 0.039367 *
M8 24974.34 4619.44 5.406 6.19e-07 ***
M9 23174.98 4615.71 5.021 2.95e-06 ***
M10 12753.36 4613.04 2.765 0.007037 **
M11 2486.49 4611.44 0.539 0.591209
t -15.63 70.15 -0.223 0.824216
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9222 on 82 degrees of freedom
Multiple R-squared: 0.802, Adjusted R-squared: 0.7706
F-statistic: 25.55 on 13 and 82 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,] 1.143745e-02 2.287491e-02 0.988562547
[2,] 4.781230e-03 9.562461e-03 0.995218770
[3,] 1.900762e-03 3.801524e-03 0.998099238
[4,] 1.356098e-03 2.712196e-03 0.998643902
[5,] 4.570893e-04 9.141786e-04 0.999542911
[6,] 1.573122e-04 3.146244e-04 0.999842688
[7,] 8.862303e-05 1.772461e-04 0.999911377
[8,] 4.364221e-05 8.728441e-05 0.999956358
[9,] 2.783452e-05 5.566904e-05 0.999972165
[10,] 1.073066e-05 2.146132e-05 0.999989269
[11,] 4.310996e-06 8.621993e-06 0.999995689
[12,] 1.524226e-06 3.048452e-06 0.999998476
[13,] 5.375756e-07 1.075151e-06 0.999999462
[14,] 4.221776e-07 8.443552e-07 0.999999578
[15,] 1.290007e-06 2.580014e-06 0.999998710
[16,] 7.888049e-06 1.577610e-05 0.999992112
[17,] 5.359456e-06 1.071891e-05 0.999994641
[18,] 6.339977e-06 1.267995e-05 0.999993660
[19,] 1.155818e-05 2.311636e-05 0.999988442
[20,] 1.561683e-05 3.123366e-05 0.999984383
[21,] 1.034419e-05 2.068837e-05 0.999989656
[22,] 6.555824e-06 1.311165e-05 0.999993444
[23,] 6.941317e-06 1.388263e-05 0.999993059
[24,] 5.769433e-06 1.153887e-05 0.999994231
[25,] 7.946968e-06 1.589394e-05 0.999992053
[26,] 5.374879e-06 1.074976e-05 0.999994625
[27,] 3.608649e-06 7.217299e-06 0.999996391
[28,] 1.947598e-06 3.895197e-06 0.999998052
[29,] 2.605599e-06 5.211198e-06 0.999997394
[30,] 4.005761e-06 8.011523e-06 0.999995994
[31,] 8.784179e-06 1.756836e-05 0.999991216
[32,] 2.465027e-05 4.930054e-05 0.999975350
[33,] 9.785937e-05 1.957187e-04 0.999902141
[34,] 3.912751e-04 7.825501e-04 0.999608725
[35,] 1.007855e-03 2.015710e-03 0.998992145
[36,] 5.039293e-03 1.007859e-02 0.994960707
[37,] 2.045180e-02 4.090359e-02 0.979548205
[38,] 5.609447e-02 1.121889e-01 0.943905526
[39,] 4.550822e-02 9.101643e-02 0.954491785
[40,] 1.285263e-01 2.570525e-01 0.871473748
[41,] 1.428826e-01 2.857653e-01 0.857117351
[42,] 1.249688e-01 2.499375e-01 0.875031236
[43,] 1.330316e-01 2.660632e-01 0.866968397
[44,] 1.233496e-01 2.466992e-01 0.876650397
[45,] 1.245325e-01 2.490650e-01 0.875467498
[46,] 1.300695e-01 2.601391e-01 0.869930464
[47,] 1.824307e-01 3.648615e-01 0.817569262
[48,] 3.331493e-01 6.662987e-01 0.666850660
[49,] 6.546033e-01 6.907934e-01 0.345396676
[50,] 7.650113e-01 4.699773e-01 0.234988673
[51,] 9.372470e-01 1.255060e-01 0.062753012
[52,] 9.364148e-01 1.271704e-01 0.063585197
[53,] 9.734505e-01 5.309892e-02 0.026549460
[54,] 9.980991e-01 3.801879e-03 0.001900939
[55,] 9.965158e-01 6.968455e-03 0.003484227
[56,] 9.957975e-01 8.405068e-03 0.004202534
[57,] 9.967749e-01 6.450188e-03 0.003225094
[58,] 9.954754e-01 9.049135e-03 0.004524567
[59,] 9.920527e-01 1.589462e-02 0.007947312
[60,] 9.849016e-01 3.019678e-02 0.015098390
[61,] 9.710744e-01 5.785124e-02 0.028925620
[62,] 9.601900e-01 7.962008e-02 0.039810039
[63,] 9.685239e-01 6.295215e-02 0.031476073
> postscript(file="/var/www/html/rcomp/tmp/1p8pr1258546290.ps",horizontal=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/2bv751258546290.ps",horizontal=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/337pj1258546290.ps",horizontal=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/4nwpi1258546290.ps",horizontal=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/5ltab1258546290.ps",horizontal=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
14618.300000 12379.925000 9291.300000 8875.800000 7105.925000
6 7 8 9 10
6654.675000 10020.800000 8784.300000 8678.300000 9659.550000
11 12 13 14 15
11604.050000 10815.175000 6727.891667 6267.516667 5657.891667
16 17 18 19 20
4957.391667 3424.516667 4650.266667 2679.391667 -275.108333
21 22 23 24 25
4599.891667 3016.141667 3671.641667 3866.766667 -597.516667
26 27 28 29 30
86.108333 -882.516667 -2483.016667 -3188.891667 -6174.141667
31 32 33 34 35
-7608.016667 -320.516667 -3058.516667 -6390.266667 -5769.766667
36 37 38 39 40
-4428.641667 -6014.925000 -6764.300000 -6185.925000 -7294.425000
41 42 43 44 45
-7233.300000 -9080.550000 -13387.425000 -10935.925000 -13199.925000
46 47 48 49 50
-15312.675000 -15330.175000 -16177.050000 -4167.575000 -4819.950000
51 52 53 54 55
-5945.575000 -4824.075000 -5239.950000 -5832.200000 -11285.075000
56 57 58 59 60
-6621.575000 -10328.575000 -9717.325000 -9909.825000 -11429.700000
61 62 63 64 65
-10795.983333 -9911.358333 -8182.983333 -6933.483333 -5261.358333
66 67 68 69 70
-5758.608333 -8161.483333 -6511.983333 -3992.983333 -1620.733333
71 72 73 74 75
-6355.233333 -6952.108333 -7002.391667 -4916.766667 -2318.391667
76 77 78 79 80
-1590.891667 6.233333 1227.983333 7046.108333 2071.608333
81 82 83 84 85
2878.608333 5664.858333 7367.358333 7669.483333 7232.200000
86 87 88 89 90
7678.825000 8566.200000 9292.700000 10386.825000 14312.575000
91 92 93 94 95
20695.700000 13809.200000 14423.200000 14700.450000 14721.950000
96
16636.075000
> postscript(file="/var/www/html/rcomp/tmp/6xwir1258546290.ps",horizontal=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 14618.300000 NA
1 12379.925000 14618.300000
2 9291.300000 12379.925000
3 8875.800000 9291.300000
4 7105.925000 8875.800000
5 6654.675000 7105.925000
6 10020.800000 6654.675000
7 8784.300000 10020.800000
8 8678.300000 8784.300000
9 9659.550000 8678.300000
10 11604.050000 9659.550000
11 10815.175000 11604.050000
12 6727.891667 10815.175000
13 6267.516667 6727.891667
14 5657.891667 6267.516667
15 4957.391667 5657.891667
16 3424.516667 4957.391667
17 4650.266667 3424.516667
18 2679.391667 4650.266667
19 -275.108333 2679.391667
20 4599.891667 -275.108333
21 3016.141667 4599.891667
22 3671.641667 3016.141667
23 3866.766667 3671.641667
24 -597.516667 3866.766667
25 86.108333 -597.516667
26 -882.516667 86.108333
27 -2483.016667 -882.516667
28 -3188.891667 -2483.016667
29 -6174.141667 -3188.891667
30 -7608.016667 -6174.141667
31 -320.516667 -7608.016667
32 -3058.516667 -320.516667
33 -6390.266667 -3058.516667
34 -5769.766667 -6390.266667
35 -4428.641667 -5769.766667
36 -6014.925000 -4428.641667
37 -6764.300000 -6014.925000
38 -6185.925000 -6764.300000
39 -7294.425000 -6185.925000
40 -7233.300000 -7294.425000
41 -9080.550000 -7233.300000
42 -13387.425000 -9080.550000
43 -10935.925000 -13387.425000
44 -13199.925000 -10935.925000
45 -15312.675000 -13199.925000
46 -15330.175000 -15312.675000
47 -16177.050000 -15330.175000
48 -4167.575000 -16177.050000
49 -4819.950000 -4167.575000
50 -5945.575000 -4819.950000
51 -4824.075000 -5945.575000
52 -5239.950000 -4824.075000
53 -5832.200000 -5239.950000
54 -11285.075000 -5832.200000
55 -6621.575000 -11285.075000
56 -10328.575000 -6621.575000
57 -9717.325000 -10328.575000
58 -9909.825000 -9717.325000
59 -11429.700000 -9909.825000
60 -10795.983333 -11429.700000
61 -9911.358333 -10795.983333
62 -8182.983333 -9911.358333
63 -6933.483333 -8182.983333
64 -5261.358333 -6933.483333
65 -5758.608333 -5261.358333
66 -8161.483333 -5758.608333
67 -6511.983333 -8161.483333
68 -3992.983333 -6511.983333
69 -1620.733333 -3992.983333
70 -6355.233333 -1620.733333
71 -6952.108333 -6355.233333
72 -7002.391667 -6952.108333
73 -4916.766667 -7002.391667
74 -2318.391667 -4916.766667
75 -1590.891667 -2318.391667
76 6.233333 -1590.891667
77 1227.983333 6.233333
78 7046.108333 1227.983333
79 2071.608333 7046.108333
80 2878.608333 2071.608333
81 5664.858333 2878.608333
82 7367.358333 5664.858333
83 7669.483333 7367.358333
84 7232.200000 7669.483333
85 7678.825000 7232.200000
86 8566.200000 7678.825000
87 9292.700000 8566.200000
88 10386.825000 9292.700000
89 14312.575000 10386.825000
90 20695.700000 14312.575000
91 13809.200000 20695.700000
92 14423.200000 13809.200000
93 14700.450000 14423.200000
94 14721.950000 14700.450000
95 16636.075000 14721.950000
96 NA 16636.075000
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 12379.925000 14618.300000
[2,] 9291.300000 12379.925000
[3,] 8875.800000 9291.300000
[4,] 7105.925000 8875.800000
[5,] 6654.675000 7105.925000
[6,] 10020.800000 6654.675000
[7,] 8784.300000 10020.800000
[8,] 8678.300000 8784.300000
[9,] 9659.550000 8678.300000
[10,] 11604.050000 9659.550000
[11,] 10815.175000 11604.050000
[12,] 6727.891667 10815.175000
[13,] 6267.516667 6727.891667
[14,] 5657.891667 6267.516667
[15,] 4957.391667 5657.891667
[16,] 3424.516667 4957.391667
[17,] 4650.266667 3424.516667
[18,] 2679.391667 4650.266667
[19,] -275.108333 2679.391667
[20,] 4599.891667 -275.108333
[21,] 3016.141667 4599.891667
[22,] 3671.641667 3016.141667
[23,] 3866.766667 3671.641667
[24,] -597.516667 3866.766667
[25,] 86.108333 -597.516667
[26,] -882.516667 86.108333
[27,] -2483.016667 -882.516667
[28,] -3188.891667 -2483.016667
[29,] -6174.141667 -3188.891667
[30,] -7608.016667 -6174.141667
[31,] -320.516667 -7608.016667
[32,] -3058.516667 -320.516667
[33,] -6390.266667 -3058.516667
[34,] -5769.766667 -6390.266667
[35,] -4428.641667 -5769.766667
[36,] -6014.925000 -4428.641667
[37,] -6764.300000 -6014.925000
[38,] -6185.925000 -6764.300000
[39,] -7294.425000 -6185.925000
[40,] -7233.300000 -7294.425000
[41,] -9080.550000 -7233.300000
[42,] -13387.425000 -9080.550000
[43,] -10935.925000 -13387.425000
[44,] -13199.925000 -10935.925000
[45,] -15312.675000 -13199.925000
[46,] -15330.175000 -15312.675000
[47,] -16177.050000 -15330.175000
[48,] -4167.575000 -16177.050000
[49,] -4819.950000 -4167.575000
[50,] -5945.575000 -4819.950000
[51,] -4824.075000 -5945.575000
[52,] -5239.950000 -4824.075000
[53,] -5832.200000 -5239.950000
[54,] -11285.075000 -5832.200000
[55,] -6621.575000 -11285.075000
[56,] -10328.575000 -6621.575000
[57,] -9717.325000 -10328.575000
[58,] -9909.825000 -9717.325000
[59,] -11429.700000 -9909.825000
[60,] -10795.983333 -11429.700000
[61,] -9911.358333 -10795.983333
[62,] -8182.983333 -9911.358333
[63,] -6933.483333 -8182.983333
[64,] -5261.358333 -6933.483333
[65,] -5758.608333 -5261.358333
[66,] -8161.483333 -5758.608333
[67,] -6511.983333 -8161.483333
[68,] -3992.983333 -6511.983333
[69,] -1620.733333 -3992.983333
[70,] -6355.233333 -1620.733333
[71,] -6952.108333 -6355.233333
[72,] -7002.391667 -6952.108333
[73,] -4916.766667 -7002.391667
[74,] -2318.391667 -4916.766667
[75,] -1590.891667 -2318.391667
[76,] 6.233333 -1590.891667
[77,] 1227.983333 6.233333
[78,] 7046.108333 1227.983333
[79,] 2071.608333 7046.108333
[80,] 2878.608333 2071.608333
[81,] 5664.858333 2878.608333
[82,] 7367.358333 5664.858333
[83,] 7669.483333 7367.358333
[84,] 7232.200000 7669.483333
[85,] 7678.825000 7232.200000
[86,] 8566.200000 7678.825000
[87,] 9292.700000 8566.200000
[88,] 10386.825000 9292.700000
[89,] 14312.575000 10386.825000
[90,] 20695.700000 14312.575000
[91,] 13809.200000 20695.700000
[92,] 14423.200000 13809.200000
[93,] 14700.450000 14423.200000
[94,] 14721.950000 14700.450000
[95,] 16636.075000 14721.950000
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 12379.925000 14618.300000
2 9291.300000 12379.925000
3 8875.800000 9291.300000
4 7105.925000 8875.800000
5 6654.675000 7105.925000
6 10020.800000 6654.675000
7 8784.300000 10020.800000
8 8678.300000 8784.300000
9 9659.550000 8678.300000
10 11604.050000 9659.550000
11 10815.175000 11604.050000
12 6727.891667 10815.175000
13 6267.516667 6727.891667
14 5657.891667 6267.516667
15 4957.391667 5657.891667
16 3424.516667 4957.391667
17 4650.266667 3424.516667
18 2679.391667 4650.266667
19 -275.108333 2679.391667
20 4599.891667 -275.108333
21 3016.141667 4599.891667
22 3671.641667 3016.141667
23 3866.766667 3671.641667
24 -597.516667 3866.766667
25 86.108333 -597.516667
26 -882.516667 86.108333
27 -2483.016667 -882.516667
28 -3188.891667 -2483.016667
29 -6174.141667 -3188.891667
30 -7608.016667 -6174.141667
31 -320.516667 -7608.016667
32 -3058.516667 -320.516667
33 -6390.266667 -3058.516667
34 -5769.766667 -6390.266667
35 -4428.641667 -5769.766667
36 -6014.925000 -4428.641667
37 -6764.300000 -6014.925000
38 -6185.925000 -6764.300000
39 -7294.425000 -6185.925000
40 -7233.300000 -7294.425000
41 -9080.550000 -7233.300000
42 -13387.425000 -9080.550000
43 -10935.925000 -13387.425000
44 -13199.925000 -10935.925000
45 -15312.675000 -13199.925000
46 -15330.175000 -15312.675000
47 -16177.050000 -15330.175000
48 -4167.575000 -16177.050000
49 -4819.950000 -4167.575000
50 -5945.575000 -4819.950000
51 -4824.075000 -5945.575000
52 -5239.950000 -4824.075000
53 -5832.200000 -5239.950000
54 -11285.075000 -5832.200000
55 -6621.575000 -11285.075000
56 -10328.575000 -6621.575000
57 -9717.325000 -10328.575000
58 -9909.825000 -9717.325000
59 -11429.700000 -9909.825000
60 -10795.983333 -11429.700000
61 -9911.358333 -10795.983333
62 -8182.983333 -9911.358333
63 -6933.483333 -8182.983333
64 -5261.358333 -6933.483333
65 -5758.608333 -5261.358333
66 -8161.483333 -5758.608333
67 -6511.983333 -8161.483333
68 -3992.983333 -6511.983333
69 -1620.733333 -3992.983333
70 -6355.233333 -1620.733333
71 -6952.108333 -6355.233333
72 -7002.391667 -6952.108333
73 -4916.766667 -7002.391667
74 -2318.391667 -4916.766667
75 -1590.891667 -2318.391667
76 6.233333 -1590.891667
77 1227.983333 6.233333
78 7046.108333 1227.983333
79 2071.608333 7046.108333
80 2878.608333 2071.608333
81 5664.858333 2878.608333
82 7367.358333 5664.858333
83 7669.483333 7367.358333
84 7232.200000 7669.483333
85 7678.825000 7232.200000
86 8566.200000 7678.825000
87 9292.700000 8566.200000
88 10386.825000 9292.700000
89 14312.575000 10386.825000
90 20695.700000 14312.575000
91 13809.200000 20695.700000
92 14423.200000 13809.200000
93 14700.450000 14423.200000
94 14721.950000 14700.450000
95 16636.075000 14721.950000
> 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/7ch001258546290.ps",horizontal=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/8wkyy1258546290.ps",horizontal=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/94rma1258546290.ps",horizontal=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/105z1v1258546290.ps",horizontal=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/11g5jm1258546290.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/12p1re1258546290.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/13t9tr1258546290.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/14tdm21258546290.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/15ari31258546290.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/169jmp1258546290.tab")
+ }
>
> system("convert tmp/1p8pr1258546290.ps tmp/1p8pr1258546290.png")
> system("convert tmp/2bv751258546290.ps tmp/2bv751258546290.png")
> system("convert tmp/337pj1258546290.ps tmp/337pj1258546290.png")
> system("convert tmp/4nwpi1258546290.ps tmp/4nwpi1258546290.png")
> system("convert tmp/5ltab1258546290.ps tmp/5ltab1258546290.png")
> system("convert tmp/6xwir1258546290.ps tmp/6xwir1258546290.png")
> system("convert tmp/7ch001258546290.ps tmp/7ch001258546290.png")
> system("convert tmp/8wkyy1258546290.ps tmp/8wkyy1258546290.png")
> system("convert tmp/94rma1258546290.ps tmp/94rma1258546290.png")
> system("convert tmp/105z1v1258546290.ps tmp/105z1v1258546290.png")
>
>
> proc.time()
user system elapsed
2.884 1.609 3.483