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(106370
+ ,100.3
+ ,109375
+ ,101.9
+ ,116476
+ ,102.1
+ ,123297
+ ,103.2
+ ,114813
+ ,103.7
+ ,117925
+ ,106.2
+ ,126466
+ ,107.7
+ ,131235
+ ,109.9
+ ,120546
+ ,111.7
+ ,123791
+ ,114.9
+ ,129813
+ ,116
+ ,133463
+ ,118.3
+ ,122987
+ ,120.4
+ ,125418
+ ,126
+ ,130199
+ ,128.1
+ ,133016
+ ,130.1
+ ,121454
+ ,130.8
+ ,122044
+ ,133.6
+ ,128313
+ ,134.2
+ ,131556
+ ,135.5
+ ,120027
+ ,136.2
+ ,123001
+ ,139.1
+ ,130111
+ ,139
+ ,132524
+ ,139.6
+ ,123742
+ ,138.7
+ ,124931
+ ,140.9
+ ,133646
+ ,141.3
+ ,136557
+ ,141.8
+ ,127509
+ ,142
+ ,128945
+ ,144.5
+ ,137191
+ ,144.6
+ ,139716
+ ,145.5
+ ,129083
+ ,146.8
+ ,131604
+ ,149.5
+ ,139413
+ ,149.9
+ ,143125
+ ,150.1
+ ,133948
+ ,150.9
+ ,137116
+ ,152.8
+ ,144864
+ ,153.1
+ ,149277
+ ,154
+ ,138796
+ ,154.9
+ ,143258
+ ,156.9
+ ,150034
+ ,158.4
+ ,154708
+ ,159.7
+ ,144888
+ ,160.2
+ ,148762
+ ,163.2
+ ,156500
+ ,163.7
+ ,161088
+ ,164.4
+ ,152772
+ ,163.7
+ ,158011
+ ,165.5
+ ,163318
+ ,165.6
+ ,169969
+ ,166.8
+ ,162269
+ ,167.5
+ ,165765
+ ,170.6
+ ,170600
+ ,170.9
+ ,174681
+ ,172
+ ,166364
+ ,171.8
+ ,170240
+ ,173.9
+ ,176150
+ ,174
+ ,182056
+ ,173.8
+ ,172218
+ ,173.9
+ ,177856
+ ,176
+ ,182253
+ ,176.6
+ ,188090
+ ,178.2
+ ,176863
+ ,179.2
+ ,183273
+ ,181.3
+ ,187969
+ ,181.8
+ ,194650
+ ,182.9
+ ,183036
+ ,183.8
+ ,189516
+ ,186.3
+ ,193805
+ ,187.4
+ ,200499
+ ,189.2
+ ,188142
+ ,189.7
+ ,193732
+ ,191.9
+ ,197126
+ ,192.6
+ ,205140
+ ,193.7
+ ,191751
+ ,194.2
+ ,196700
+ ,197.6
+ ,199784
+ ,199.3
+ ,207360
+ ,201.4
+ ,196101
+ ,203
+ ,200824
+ ,206.3
+ ,205743
+ ,207.1
+ ,212489
+ ,209.8
+ ,200810
+ ,211.1
+ ,203683
+ ,215.3
+ ,207286
+ ,217.4
+ ,210910
+ ,215.5
+ ,194915
+ ,210.9
+ ,217920
+ ,212.6)
+ ,dim=c(2
+ ,90)
+ ,dimnames=list(c('Y'
+ ,'X')
+ ,1:90))
> y <- array(NA,dim=c(2,90),dimnames=list(c('Y','X'),1:90))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 106370 100.3 1 0 0 0 0 0 0 0 0 0 0
2 109375 101.9 0 1 0 0 0 0 0 0 0 0 0
3 116476 102.1 0 0 1 0 0 0 0 0 0 0 0
4 123297 103.2 0 0 0 1 0 0 0 0 0 0 0
5 114813 103.7 0 0 0 0 1 0 0 0 0 0 0
6 117925 106.2 0 0 0 0 0 1 0 0 0 0 0
7 126466 107.7 0 0 0 0 0 0 1 0 0 0 0
8 131235 109.9 0 0 0 0 0 0 0 1 0 0 0
9 120546 111.7 0 0 0 0 0 0 0 0 1 0 0
10 123791 114.9 0 0 0 0 0 0 0 0 0 1 0
11 129813 116.0 0 0 0 0 0 0 0 0 0 0 1
12 133463 118.3 0 0 0 0 0 0 0 0 0 0 0
13 122987 120.4 1 0 0 0 0 0 0 0 0 0 0
14 125418 126.0 0 1 0 0 0 0 0 0 0 0 0
15 130199 128.1 0 0 1 0 0 0 0 0 0 0 0
16 133016 130.1 0 0 0 1 0 0 0 0 0 0 0
17 121454 130.8 0 0 0 0 1 0 0 0 0 0 0
18 122044 133.6 0 0 0 0 0 1 0 0 0 0 0
19 128313 134.2 0 0 0 0 0 0 1 0 0 0 0
20 131556 135.5 0 0 0 0 0 0 0 1 0 0 0
21 120027 136.2 0 0 0 0 0 0 0 0 1 0 0
22 123001 139.1 0 0 0 0 0 0 0 0 0 1 0
23 130111 139.0 0 0 0 0 0 0 0 0 0 0 1
24 132524 139.6 0 0 0 0 0 0 0 0 0 0 0
25 123742 138.7 1 0 0 0 0 0 0 0 0 0 0
26 124931 140.9 0 1 0 0 0 0 0 0 0 0 0
27 133646 141.3 0 0 1 0 0 0 0 0 0 0 0
28 136557 141.8 0 0 0 1 0 0 0 0 0 0 0
29 127509 142.0 0 0 0 0 1 0 0 0 0 0 0
30 128945 144.5 0 0 0 0 0 1 0 0 0 0 0
31 137191 144.6 0 0 0 0 0 0 1 0 0 0 0
32 139716 145.5 0 0 0 0 0 0 0 1 0 0 0
33 129083 146.8 0 0 0 0 0 0 0 0 1 0 0
34 131604 149.5 0 0 0 0 0 0 0 0 0 1 0
35 139413 149.9 0 0 0 0 0 0 0 0 0 0 1
36 143125 150.1 0 0 0 0 0 0 0 0 0 0 0
37 133948 150.9 1 0 0 0 0 0 0 0 0 0 0
38 137116 152.8 0 1 0 0 0 0 0 0 0 0 0
39 144864 153.1 0 0 1 0 0 0 0 0 0 0 0
40 149277 154.0 0 0 0 1 0 0 0 0 0 0 0
41 138796 154.9 0 0 0 0 1 0 0 0 0 0 0
42 143258 156.9 0 0 0 0 0 1 0 0 0 0 0
43 150034 158.4 0 0 0 0 0 0 1 0 0 0 0
44 154708 159.7 0 0 0 0 0 0 0 1 0 0 0
45 144888 160.2 0 0 0 0 0 0 0 0 1 0 0
46 148762 163.2 0 0 0 0 0 0 0 0 0 1 0
47 156500 163.7 0 0 0 0 0 0 0 0 0 0 1
48 161088 164.4 0 0 0 0 0 0 0 0 0 0 0
49 152772 163.7 1 0 0 0 0 0 0 0 0 0 0
50 158011 165.5 0 1 0 0 0 0 0 0 0 0 0
51 163318 165.6 0 0 1 0 0 0 0 0 0 0 0
52 169969 166.8 0 0 0 1 0 0 0 0 0 0 0
53 162269 167.5 0 0 0 0 1 0 0 0 0 0 0
54 165765 170.6 0 0 0 0 0 1 0 0 0 0 0
55 170600 170.9 0 0 0 0 0 0 1 0 0 0 0
56 174681 172.0 0 0 0 0 0 0 0 1 0 0 0
57 166364 171.8 0 0 0 0 0 0 0 0 1 0 0
58 170240 173.9 0 0 0 0 0 0 0 0 0 1 0
59 176150 174.0 0 0 0 0 0 0 0 0 0 0 1
60 182056 173.8 0 0 0 0 0 0 0 0 0 0 0
61 172218 173.9 1 0 0 0 0 0 0 0 0 0 0
62 177856 176.0 0 1 0 0 0 0 0 0 0 0 0
63 182253 176.6 0 0 1 0 0 0 0 0 0 0 0
64 188090 178.2 0 0 0 1 0 0 0 0 0 0 0
65 176863 179.2 0 0 0 0 1 0 0 0 0 0 0
66 183273 181.3 0 0 0 0 0 1 0 0 0 0 0
67 187969 181.8 0 0 0 0 0 0 1 0 0 0 0
68 194650 182.9 0 0 0 0 0 0 0 1 0 0 0
69 183036 183.8 0 0 0 0 0 0 0 0 1 0 0
70 189516 186.3 0 0 0 0 0 0 0 0 0 1 0
71 193805 187.4 0 0 0 0 0 0 0 0 0 0 1
72 200499 189.2 0 0 0 0 0 0 0 0 0 0 0
73 188142 189.7 1 0 0 0 0 0 0 0 0 0 0
74 193732 191.9 0 1 0 0 0 0 0 0 0 0 0
75 197126 192.6 0 0 1 0 0 0 0 0 0 0 0
76 205140 193.7 0 0 0 1 0 0 0 0 0 0 0
77 191751 194.2 0 0 0 0 1 0 0 0 0 0 0
78 196700 197.6 0 0 0 0 0 1 0 0 0 0 0
79 199784 199.3 0 0 0 0 0 0 1 0 0 0 0
80 207360 201.4 0 0 0 0 0 0 0 1 0 0 0
81 196101 203.0 0 0 0 0 0 0 0 0 1 0 0
82 200824 206.3 0 0 0 0 0 0 0 0 0 1 0
83 205743 207.1 0 0 0 0 0 0 0 0 0 0 1
84 212489 209.8 0 0 0 0 0 0 0 0 0 0 0
85 200810 211.1 1 0 0 0 0 0 0 0 0 0 0
86 203683 215.3 0 1 0 0 0 0 0 0 0 0 0
87 207286 217.4 0 0 1 0 0 0 0 0 0 0 0
88 210910 215.5 0 0 0 1 0 0 0 0 0 0 0
89 194915 210.9 0 0 0 0 1 0 0 0 0 0 0
90 217920 212.6 0 0 0 0 0 1 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
12578.3 940.6 -9273.4 -8171.5 -3305.0 1066.8
M5 M6 M7 M8 M9 M10
-9907.2 -6338.0 -2779.3 669.7 -10768.8 -9459.8
M11
-3727.1
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-13214.6 -8702.8 240.9 6856.6 15362.4
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12578.35 6086.45 2.067 0.0421 *
X 940.62 30.69 30.654 <2e-16 ***
M1 -9273.41 4718.07 -1.966 0.0530 .
M2 -8171.46 4714.75 -1.733 0.0871 .
M3 -3304.96 4714.03 -0.701 0.4854
M4 1066.79 4713.45 0.226 0.8215
M5 -9907.20 4713.46 -2.102 0.0388 *
M6 -6338.00 4712.48 -1.345 0.1826
M7 -2779.31 4871.58 -0.571 0.5700
M8 669.67 4869.87 0.138 0.8910
M9 -10768.77 4868.96 -2.212 0.0300 *
M10 -9459.80 4867.26 -1.944 0.0556 .
M11 -3727.14 4867.11 -0.766 0.4461
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9105 on 77 degrees of freedom
Multiple R-squared: 0.9262, Adjusted R-squared: 0.9147
F-statistic: 80.55 on 12 and 77 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.179096517 0.3581930335 0.8209034833
[2,] 0.214837660 0.4296753195 0.7851623403
[3,] 0.249279036 0.4985580712 0.7507209644
[4,] 0.276457763 0.5529155266 0.7235422367
[5,] 0.288868523 0.5777370468 0.7111314766
[6,] 0.282237804 0.5644756079 0.7177621960
[7,] 0.262387504 0.5247750076 0.7376124962
[8,] 0.213807801 0.4276156020 0.7861921990
[9,] 0.179303353 0.3586067055 0.8206966472
[10,] 0.127845630 0.2556912591 0.8721543704
[11,] 0.086065658 0.1721313169 0.9139343416
[12,] 0.062928156 0.1258563114 0.9370718443
[13,] 0.040795878 0.0815917564 0.9592041218
[14,] 0.026818170 0.0536363401 0.9731818299
[15,] 0.018572938 0.0371458768 0.9814270616
[16,] 0.012260530 0.0245210604 0.9877394698
[17,] 0.007598492 0.0151969847 0.9924015076
[18,] 0.004895161 0.0097903226 0.9951048387
[19,] 0.003344450 0.0066889000 0.9966555500
[20,] 0.002268201 0.0045364021 0.9977317990
[21,] 0.001805173 0.0036103464 0.9981948268
[22,] 0.002325351 0.0046507011 0.9976746495
[23,] 0.003695259 0.0073905186 0.9963047407
[24,] 0.005199146 0.0103982925 0.9948008538
[25,] 0.007418518 0.0148370370 0.9925814815
[26,] 0.008543664 0.0170873271 0.9914563365
[27,] 0.020356605 0.0407132095 0.9796433952
[28,] 0.028696966 0.0573939324 0.9713030338
[29,] 0.049465136 0.0989302723 0.9505348639
[30,] 0.088298426 0.1765968515 0.9117015742
[31,] 0.178567858 0.3571357151 0.8214321424
[32,] 0.282814914 0.5656298280 0.7171850860
[33,] 0.474005198 0.9480103966 0.5259948017
[34,] 0.639751977 0.7204960464 0.3602480232
[35,] 0.794523233 0.4109535332 0.2054767666
[36,] 0.852369292 0.2952614167 0.1476307083
[37,] 0.907061568 0.1858768632 0.0929384316
[38,] 0.930884648 0.1382307032 0.0691153516
[39,] 0.979582983 0.0408340331 0.0204170165
[40,] 0.986478122 0.0270437552 0.0135218776
[41,] 0.994988462 0.0100230757 0.0050115379
[42,] 0.996972332 0.0060553366 0.0030276683
[43,] 0.998720291 0.0025594175 0.0012797087
[44,] 0.999198636 0.0016027285 0.0008013642
[45,] 0.999524081 0.0009518374 0.0004759187
[46,] 0.999493703 0.0010125943 0.0005062971
[47,] 0.999350811 0.0012983788 0.0006491894
[48,] 0.998902960 0.0021940808 0.0010970404
[49,] 0.998332573 0.0033348544 0.0016674272
[50,] 0.996905866 0.0061882682 0.0030941341
[51,] 0.999066307 0.0018673857 0.0009336929
[52,] 0.998020081 0.0039598382 0.0019799191
[53,] 0.996153980 0.0076920397 0.0038460199
[54,] 0.992271357 0.0154572856 0.0077286428
[55,] 0.983638558 0.0327228842 0.0163614421
[56,] 0.965281281 0.0694374381 0.0347187190
[57,] 0.929051974 0.1418960524 0.0709480262
[58,] 0.859148523 0.2817029531 0.1408514766
[59,] 0.730610403 0.5387791931 0.2693895965
> postscript(file="/var/www/html/rcomp/tmp/1p48j1258926973.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/2sgf71258926973.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/3biau1258926973.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/45vie1258926973.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/5eqv31258926973.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 = 90
Frequency = 1
1 2 3 4 5 6
8721.09385 9119.14850 11165.52694 12580.09930 14599.78265 11791.04037
7 8 9 10 11 12
15362.41799 14613.08421 13669.41176 12595.45915 11850.12374 9609.56048
13 14 15 16 17 18
6431.67487 2493.25808 432.46258 -3003.52113 -4249.96135 -9862.88899
19 20 21 22 23 24
-7716.95530 -9145.73300 -9894.72581 -10957.49305 -9486.08703 -11364.59994
25 26 27 28 29 30
-10026.63196 -12008.94803 -8536.69317 -10467.75009 -8729.88138 -13214.62366
31 32 33 34 35 36
-8621.38104 -10391.91160 -10809.27512 -12136.91879 -10436.82170 -10640.08746
37 38 39 40 41 42
-11296.16985 -11017.30056 -8417.98391 -9223.28797 -9576.85177 -10565.28512
43 44 45 46 47 48
-8758.90750 -8756.68521 -7608.55444 -7865.38347 -6330.34817 -6127.92286
49 50 51 52 53 54
-4512.07846 -2068.14738 -1721.70716 -571.19658 2044.36320 -944.74980
55 56 57 58 59 60
49.36925 -353.28488 2956.27838 3548.00543 3631.28788 5998.26926
61 62 63 64 65 66
5339.61937 7900.36509 6866.49638 6826.75982 5633.13424 6498.63910
67 68 69 70 71 72
7165.63458 9362.98045 8340.86407 11160.34397 8682.00856 9955.75422
73 74 75 76 77 78
6401.85719 8820.54112 6689.61063 9297.18299 6411.86634 4593.56799
79 80 81 82 83 84
2519.82203 4671.55004 3346.00116 3655.98677 2089.83672 2569.02631
85 86 87 88 89 90
-1059.36501 -3238.91680 -6477.71229 -5438.28635 -6132.45192 11704.30009
> postscript(file="/var/www/html/rcomp/tmp/6pqpc1258926973.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 = 90
Frequency = 1
lag(myerror, k = 1) myerror
0 8721.09385 NA
1 9119.14850 8721.09385
2 11165.52694 9119.14850
3 12580.09930 11165.52694
4 14599.78265 12580.09930
5 11791.04037 14599.78265
6 15362.41799 11791.04037
7 14613.08421 15362.41799
8 13669.41176 14613.08421
9 12595.45915 13669.41176
10 11850.12374 12595.45915
11 9609.56048 11850.12374
12 6431.67487 9609.56048
13 2493.25808 6431.67487
14 432.46258 2493.25808
15 -3003.52113 432.46258
16 -4249.96135 -3003.52113
17 -9862.88899 -4249.96135
18 -7716.95530 -9862.88899
19 -9145.73300 -7716.95530
20 -9894.72581 -9145.73300
21 -10957.49305 -9894.72581
22 -9486.08703 -10957.49305
23 -11364.59994 -9486.08703
24 -10026.63196 -11364.59994
25 -12008.94803 -10026.63196
26 -8536.69317 -12008.94803
27 -10467.75009 -8536.69317
28 -8729.88138 -10467.75009
29 -13214.62366 -8729.88138
30 -8621.38104 -13214.62366
31 -10391.91160 -8621.38104
32 -10809.27512 -10391.91160
33 -12136.91879 -10809.27512
34 -10436.82170 -12136.91879
35 -10640.08746 -10436.82170
36 -11296.16985 -10640.08746
37 -11017.30056 -11296.16985
38 -8417.98391 -11017.30056
39 -9223.28797 -8417.98391
40 -9576.85177 -9223.28797
41 -10565.28512 -9576.85177
42 -8758.90750 -10565.28512
43 -8756.68521 -8758.90750
44 -7608.55444 -8756.68521
45 -7865.38347 -7608.55444
46 -6330.34817 -7865.38347
47 -6127.92286 -6330.34817
48 -4512.07846 -6127.92286
49 -2068.14738 -4512.07846
50 -1721.70716 -2068.14738
51 -571.19658 -1721.70716
52 2044.36320 -571.19658
53 -944.74980 2044.36320
54 49.36925 -944.74980
55 -353.28488 49.36925
56 2956.27838 -353.28488
57 3548.00543 2956.27838
58 3631.28788 3548.00543
59 5998.26926 3631.28788
60 5339.61937 5998.26926
61 7900.36509 5339.61937
62 6866.49638 7900.36509
63 6826.75982 6866.49638
64 5633.13424 6826.75982
65 6498.63910 5633.13424
66 7165.63458 6498.63910
67 9362.98045 7165.63458
68 8340.86407 9362.98045
69 11160.34397 8340.86407
70 8682.00856 11160.34397
71 9955.75422 8682.00856
72 6401.85719 9955.75422
73 8820.54112 6401.85719
74 6689.61063 8820.54112
75 9297.18299 6689.61063
76 6411.86634 9297.18299
77 4593.56799 6411.86634
78 2519.82203 4593.56799
79 4671.55004 2519.82203
80 3346.00116 4671.55004
81 3655.98677 3346.00116
82 2089.83672 3655.98677
83 2569.02631 2089.83672
84 -1059.36501 2569.02631
85 -3238.91680 -1059.36501
86 -6477.71229 -3238.91680
87 -5438.28635 -6477.71229
88 -6132.45192 -5438.28635
89 11704.30009 -6132.45192
90 NA 11704.30009
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 9119.14850 8721.09385
[2,] 11165.52694 9119.14850
[3,] 12580.09930 11165.52694
[4,] 14599.78265 12580.09930
[5,] 11791.04037 14599.78265
[6,] 15362.41799 11791.04037
[7,] 14613.08421 15362.41799
[8,] 13669.41176 14613.08421
[9,] 12595.45915 13669.41176
[10,] 11850.12374 12595.45915
[11,] 9609.56048 11850.12374
[12,] 6431.67487 9609.56048
[13,] 2493.25808 6431.67487
[14,] 432.46258 2493.25808
[15,] -3003.52113 432.46258
[16,] -4249.96135 -3003.52113
[17,] -9862.88899 -4249.96135
[18,] -7716.95530 -9862.88899
[19,] -9145.73300 -7716.95530
[20,] -9894.72581 -9145.73300
[21,] -10957.49305 -9894.72581
[22,] -9486.08703 -10957.49305
[23,] -11364.59994 -9486.08703
[24,] -10026.63196 -11364.59994
[25,] -12008.94803 -10026.63196
[26,] -8536.69317 -12008.94803
[27,] -10467.75009 -8536.69317
[28,] -8729.88138 -10467.75009
[29,] -13214.62366 -8729.88138
[30,] -8621.38104 -13214.62366
[31,] -10391.91160 -8621.38104
[32,] -10809.27512 -10391.91160
[33,] -12136.91879 -10809.27512
[34,] -10436.82170 -12136.91879
[35,] -10640.08746 -10436.82170
[36,] -11296.16985 -10640.08746
[37,] -11017.30056 -11296.16985
[38,] -8417.98391 -11017.30056
[39,] -9223.28797 -8417.98391
[40,] -9576.85177 -9223.28797
[41,] -10565.28512 -9576.85177
[42,] -8758.90750 -10565.28512
[43,] -8756.68521 -8758.90750
[44,] -7608.55444 -8756.68521
[45,] -7865.38347 -7608.55444
[46,] -6330.34817 -7865.38347
[47,] -6127.92286 -6330.34817
[48,] -4512.07846 -6127.92286
[49,] -2068.14738 -4512.07846
[50,] -1721.70716 -2068.14738
[51,] -571.19658 -1721.70716
[52,] 2044.36320 -571.19658
[53,] -944.74980 2044.36320
[54,] 49.36925 -944.74980
[55,] -353.28488 49.36925
[56,] 2956.27838 -353.28488
[57,] 3548.00543 2956.27838
[58,] 3631.28788 3548.00543
[59,] 5998.26926 3631.28788
[60,] 5339.61937 5998.26926
[61,] 7900.36509 5339.61937
[62,] 6866.49638 7900.36509
[63,] 6826.75982 6866.49638
[64,] 5633.13424 6826.75982
[65,] 6498.63910 5633.13424
[66,] 7165.63458 6498.63910
[67,] 9362.98045 7165.63458
[68,] 8340.86407 9362.98045
[69,] 11160.34397 8340.86407
[70,] 8682.00856 11160.34397
[71,] 9955.75422 8682.00856
[72,] 6401.85719 9955.75422
[73,] 8820.54112 6401.85719
[74,] 6689.61063 8820.54112
[75,] 9297.18299 6689.61063
[76,] 6411.86634 9297.18299
[77,] 4593.56799 6411.86634
[78,] 2519.82203 4593.56799
[79,] 4671.55004 2519.82203
[80,] 3346.00116 4671.55004
[81,] 3655.98677 3346.00116
[82,] 2089.83672 3655.98677
[83,] 2569.02631 2089.83672
[84,] -1059.36501 2569.02631
[85,] -3238.91680 -1059.36501
[86,] -6477.71229 -3238.91680
[87,] -5438.28635 -6477.71229
[88,] -6132.45192 -5438.28635
[89,] 11704.30009 -6132.45192
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 9119.14850 8721.09385
2 11165.52694 9119.14850
3 12580.09930 11165.52694
4 14599.78265 12580.09930
5 11791.04037 14599.78265
6 15362.41799 11791.04037
7 14613.08421 15362.41799
8 13669.41176 14613.08421
9 12595.45915 13669.41176
10 11850.12374 12595.45915
11 9609.56048 11850.12374
12 6431.67487 9609.56048
13 2493.25808 6431.67487
14 432.46258 2493.25808
15 -3003.52113 432.46258
16 -4249.96135 -3003.52113
17 -9862.88899 -4249.96135
18 -7716.95530 -9862.88899
19 -9145.73300 -7716.95530
20 -9894.72581 -9145.73300
21 -10957.49305 -9894.72581
22 -9486.08703 -10957.49305
23 -11364.59994 -9486.08703
24 -10026.63196 -11364.59994
25 -12008.94803 -10026.63196
26 -8536.69317 -12008.94803
27 -10467.75009 -8536.69317
28 -8729.88138 -10467.75009
29 -13214.62366 -8729.88138
30 -8621.38104 -13214.62366
31 -10391.91160 -8621.38104
32 -10809.27512 -10391.91160
33 -12136.91879 -10809.27512
34 -10436.82170 -12136.91879
35 -10640.08746 -10436.82170
36 -11296.16985 -10640.08746
37 -11017.30056 -11296.16985
38 -8417.98391 -11017.30056
39 -9223.28797 -8417.98391
40 -9576.85177 -9223.28797
41 -10565.28512 -9576.85177
42 -8758.90750 -10565.28512
43 -8756.68521 -8758.90750
44 -7608.55444 -8756.68521
45 -7865.38347 -7608.55444
46 -6330.34817 -7865.38347
47 -6127.92286 -6330.34817
48 -4512.07846 -6127.92286
49 -2068.14738 -4512.07846
50 -1721.70716 -2068.14738
51 -571.19658 -1721.70716
52 2044.36320 -571.19658
53 -944.74980 2044.36320
54 49.36925 -944.74980
55 -353.28488 49.36925
56 2956.27838 -353.28488
57 3548.00543 2956.27838
58 3631.28788 3548.00543
59 5998.26926 3631.28788
60 5339.61937 5998.26926
61 7900.36509 5339.61937
62 6866.49638 7900.36509
63 6826.75982 6866.49638
64 5633.13424 6826.75982
65 6498.63910 5633.13424
66 7165.63458 6498.63910
67 9362.98045 7165.63458
68 8340.86407 9362.98045
69 11160.34397 8340.86407
70 8682.00856 11160.34397
71 9955.75422 8682.00856
72 6401.85719 9955.75422
73 8820.54112 6401.85719
74 6689.61063 8820.54112
75 9297.18299 6689.61063
76 6411.86634 9297.18299
77 4593.56799 6411.86634
78 2519.82203 4593.56799
79 4671.55004 2519.82203
80 3346.00116 4671.55004
81 3655.98677 3346.00116
82 2089.83672 3655.98677
83 2569.02631 2089.83672
84 -1059.36501 2569.02631
85 -3238.91680 -1059.36501
86 -6477.71229 -3238.91680
87 -5438.28635 -6477.71229
88 -6132.45192 -5438.28635
89 11704.30009 -6132.45192
> 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/78pt51258926973.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/8d6y01258926973.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/9mkb01258926973.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/10rtlc1258926973.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/112ivw1258926973.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/12cxsf1258926973.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/13fdmi1258926973.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/14brsj1258926973.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/15udba1258926973.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/1673r31258926973.tab")
+ }
>
> system("convert tmp/1p48j1258926973.ps tmp/1p48j1258926973.png")
> system("convert tmp/2sgf71258926973.ps tmp/2sgf71258926973.png")
> system("convert tmp/3biau1258926973.ps tmp/3biau1258926973.png")
> system("convert tmp/45vie1258926973.ps tmp/45vie1258926973.png")
> system("convert tmp/5eqv31258926973.ps tmp/5eqv31258926973.png")
> system("convert tmp/6pqpc1258926973.ps tmp/6pqpc1258926973.png")
> system("convert tmp/78pt51258926973.ps tmp/78pt51258926973.png")
> system("convert tmp/8d6y01258926973.ps tmp/8d6y01258926973.png")
> system("convert tmp/9mkb01258926973.ps tmp/9mkb01258926973.png")
> system("convert tmp/10rtlc1258926973.ps tmp/10rtlc1258926973.png")
>
>
> proc.time()
user system elapsed
2.865 1.601 4.047