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(274412
+ ,244752
+ ,272433
+ ,244576
+ ,268361
+ ,241572
+ ,268586
+ ,240541
+ ,264768
+ ,236089
+ ,269974
+ ,236997
+ ,304744
+ ,264579
+ ,309365
+ ,270349
+ ,308347
+ ,269645
+ ,298427
+ ,267037
+ ,289231
+ ,258113
+ ,291975
+ ,262813
+ ,294912
+ ,267413
+ ,293488
+ ,267366
+ ,290555
+ ,264777
+ ,284736
+ ,258863
+ ,281818
+ ,254844
+ ,287854
+ ,254868
+ ,316263
+ ,277267
+ ,325412
+ ,285351
+ ,326011
+ ,286602
+ ,328282
+ ,283042
+ ,317480
+ ,276687
+ ,317539
+ ,277915
+ ,313737
+ ,277128
+ ,312276
+ ,277103
+ ,309391
+ ,275037
+ ,302950
+ ,270150
+ ,300316
+ ,267140
+ ,304035
+ ,264993
+ ,333476
+ ,287259
+ ,337698
+ ,291186
+ ,335932
+ ,292300
+ ,323931
+ ,288186
+ ,313927
+ ,281477
+ ,314485
+ ,282656
+ ,313218
+ ,280190
+ ,309664
+ ,280408
+ ,302963
+ ,276836
+ ,298989
+ ,275216
+ ,298423
+ ,274352
+ ,301631
+ ,271311
+ ,329765
+ ,289802
+ ,335083
+ ,290726
+ ,327616
+ ,292300
+ ,309119
+ ,278506
+ ,295916
+ ,269826
+ ,291413
+ ,265861
+ ,291542
+ ,269034
+ ,284678
+ ,264176
+ ,276475
+ ,255198
+ ,272566
+ ,253353
+ ,264981
+ ,246057
+ ,263290
+ ,235372
+ ,296806
+ ,258556
+ ,303598
+ ,260993
+ ,286994
+ ,254663
+ ,276427
+ ,250643
+ ,266424
+ ,243422
+ ,267153
+ ,247105
+ ,268381
+ ,248541
+ ,262522
+ ,245039
+ ,255542
+ ,237080
+ ,253158
+ ,237085
+ ,243803
+ ,225554
+ ,250741
+ ,226839
+ ,280445
+ ,247934
+ ,285257
+ ,248333
+ ,270976
+ ,246969
+ ,261076
+ ,245098
+ ,255603
+ ,246263
+ ,260376
+ ,255765
+ ,263903
+ ,264319
+ ,264291
+ ,268347
+ ,263276
+ ,273046
+ ,262572
+ ,273963
+ ,256167
+ ,267430
+ ,264221
+ ,271993
+ ,293860
+ ,292710)
+ ,dim=c(2
+ ,79)
+ ,dimnames=list(c('Y'
+ ,'X')
+ ,1:79))
> y <- array(NA,dim=c(2,79),dimnames=list(c('Y','X'),1:79))
> 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 274412 244752 1 0 0 0 0 0 0 0 0 0 0
2 272433 244576 0 1 0 0 0 0 0 0 0 0 0
3 268361 241572 0 0 1 0 0 0 0 0 0 0 0
4 268586 240541 0 0 0 1 0 0 0 0 0 0 0
5 264768 236089 0 0 0 0 1 0 0 0 0 0 0
6 269974 236997 0 0 0 0 0 1 0 0 0 0 0
7 304744 264579 0 0 0 0 0 0 1 0 0 0 0
8 309365 270349 0 0 0 0 0 0 0 1 0 0 0
9 308347 269645 0 0 0 0 0 0 0 0 1 0 0
10 298427 267037 0 0 0 0 0 0 0 0 0 1 0
11 289231 258113 0 0 0 0 0 0 0 0 0 0 1
12 291975 262813 0 0 0 0 0 0 0 0 0 0 0
13 294912 267413 1 0 0 0 0 0 0 0 0 0 0
14 293488 267366 0 1 0 0 0 0 0 0 0 0 0
15 290555 264777 0 0 1 0 0 0 0 0 0 0 0
16 284736 258863 0 0 0 1 0 0 0 0 0 0 0
17 281818 254844 0 0 0 0 1 0 0 0 0 0 0
18 287854 254868 0 0 0 0 0 1 0 0 0 0 0
19 316263 277267 0 0 0 0 0 0 1 0 0 0 0
20 325412 285351 0 0 0 0 0 0 0 1 0 0 0
21 326011 286602 0 0 0 0 0 0 0 0 1 0 0
22 328282 283042 0 0 0 0 0 0 0 0 0 1 0
23 317480 276687 0 0 0 0 0 0 0 0 0 0 1
24 317539 277915 0 0 0 0 0 0 0 0 0 0 0
25 313737 277128 1 0 0 0 0 0 0 0 0 0 0
26 312276 277103 0 1 0 0 0 0 0 0 0 0 0
27 309391 275037 0 0 1 0 0 0 0 0 0 0 0
28 302950 270150 0 0 0 1 0 0 0 0 0 0 0
29 300316 267140 0 0 0 0 1 0 0 0 0 0 0
30 304035 264993 0 0 0 0 0 1 0 0 0 0 0
31 333476 287259 0 0 0 0 0 0 1 0 0 0 0
32 337698 291186 0 0 0 0 0 0 0 1 0 0 0
33 335932 292300 0 0 0 0 0 0 0 0 1 0 0
34 323931 288186 0 0 0 0 0 0 0 0 0 1 0
35 313927 281477 0 0 0 0 0 0 0 0 0 0 1
36 314485 282656 0 0 0 0 0 0 0 0 0 0 0
37 313218 280190 1 0 0 0 0 0 0 0 0 0 0
38 309664 280408 0 1 0 0 0 0 0 0 0 0 0
39 302963 276836 0 0 1 0 0 0 0 0 0 0 0
40 298989 275216 0 0 0 1 0 0 0 0 0 0 0
41 298423 274352 0 0 0 0 1 0 0 0 0 0 0
42 301631 271311 0 0 0 0 0 1 0 0 0 0 0
43 329765 289802 0 0 0 0 0 0 1 0 0 0 0
44 335083 290726 0 0 0 0 0 0 0 1 0 0 0
45 327616 292300 0 0 0 0 0 0 0 0 1 0 0
46 309119 278506 0 0 0 0 0 0 0 0 0 1 0
47 295916 269826 0 0 0 0 0 0 0 0 0 0 1
48 291413 265861 0 0 0 0 0 0 0 0 0 0 0
49 291542 269034 1 0 0 0 0 0 0 0 0 0 0
50 284678 264176 0 1 0 0 0 0 0 0 0 0 0
51 276475 255198 0 0 1 0 0 0 0 0 0 0 0
52 272566 253353 0 0 0 1 0 0 0 0 0 0 0
53 264981 246057 0 0 0 0 1 0 0 0 0 0 0
54 263290 235372 0 0 0 0 0 1 0 0 0 0 0
55 296806 258556 0 0 0 0 0 0 1 0 0 0 0
56 303598 260993 0 0 0 0 0 0 0 1 0 0 0
57 286994 254663 0 0 0 0 0 0 0 0 1 0 0
58 276427 250643 0 0 0 0 0 0 0 0 0 1 0
59 266424 243422 0 0 0 0 0 0 0 0 0 0 1
60 267153 247105 0 0 0 0 0 0 0 0 0 0 0
61 268381 248541 1 0 0 0 0 0 0 0 0 0 0
62 262522 245039 0 1 0 0 0 0 0 0 0 0 0
63 255542 237080 0 0 1 0 0 0 0 0 0 0 0
64 253158 237085 0 0 0 1 0 0 0 0 0 0 0
65 243803 225554 0 0 0 0 1 0 0 0 0 0 0
66 250741 226839 0 0 0 0 0 1 0 0 0 0 0
67 280445 247934 0 0 0 0 0 0 1 0 0 0 0
68 285257 248333 0 0 0 0 0 0 0 1 0 0 0
69 270976 246969 0 0 0 0 0 0 0 0 1 0 0
70 261076 245098 0 0 0 0 0 0 0 0 0 1 0
71 255603 246263 0 0 0 0 0 0 0 0 0 0 1
72 260376 255765 0 0 0 0 0 0 0 0 0 0 0
73 263903 264319 1 0 0 0 0 0 0 0 0 0 0
74 264291 268347 0 1 0 0 0 0 0 0 0 0 0
75 263276 273046 0 0 1 0 0 0 0 0 0 0 0
76 262572 273963 0 0 0 1 0 0 0 0 0 0 0
77 256167 267430 0 0 0 0 1 0 0 0 0 0 0
78 264221 271993 0 0 0 0 0 1 0 0 0 0 0
79 293860 292710 0 0 0 0 0 0 1 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
3078.327 1.083 -961.336 -3251.104 -4303.814 -5366.096
M5 M6 M7 M8 M9 M10
-4286.310 1616.393 8035.320 15681.907 9730.688 5371.392
M11
2220.717
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-35078 -1883 3071 6819 13442
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.078e+03 2.506e+04 0.123 0.9026
X 1.083e+00 9.247e-02 11.714 <2e-16 ***
M1 -9.613e+02 6.976e+03 -0.138 0.8908
M2 -3.251e+03 6.976e+03 -0.466 0.6427
M3 -4.304e+03 6.989e+03 -0.616 0.5402
M4 -5.366e+03 7.004e+03 -0.766 0.4463
M5 -4.286e+03 7.067e+03 -0.607 0.5462
M6 1.616e+03 7.087e+03 0.228 0.8203
M7 8.035e+03 7.021e+03 1.144 0.2566
M8 1.568e+04 7.288e+03 2.152 0.0351 *
M9 9.731e+03 7.280e+03 1.337 0.1859
M10 5.371e+03 7.245e+03 0.741 0.4611
M11 2.221e+03 7.243e+03 0.307 0.7601
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 12540 on 66 degrees of freedom
Multiple R-squared: 0.7743, Adjusted R-squared: 0.7332
F-statistic: 18.87 on 12 and 66 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,] 5.720713e-05 1.144143e-04 0.999942793
[2,] 1.703426e-06 3.406852e-06 0.999998297
[3,] 3.745856e-07 7.491712e-07 0.999999625
[4,] 1.428398e-08 2.856795e-08 0.999999986
[5,] 8.256094e-09 1.651219e-08 0.999999992
[6,] 1.353405e-09 2.706809e-09 0.999999999
[7,] 4.632936e-05 9.265872e-05 0.999953671
[8,] 5.107729e-05 1.021546e-04 0.999948923
[9,] 5.332741e-05 1.066548e-04 0.999946673
[10,] 2.921122e-05 5.842245e-05 0.999970789
[11,] 1.517490e-05 3.034981e-05 0.999984825
[12,] 7.277383e-06 1.455477e-05 0.999992723
[13,] 2.838565e-06 5.677131e-06 0.999997161
[14,] 1.059915e-06 2.119829e-06 0.999998940
[15,] 4.346529e-07 8.693059e-07 0.999999565
[16,] 1.910325e-07 3.820650e-07 0.999999809
[17,] 7.204059e-08 1.440812e-07 0.999999928
[18,] 2.252102e-08 4.504205e-08 0.999999977
[19,] 1.286230e-08 2.572459e-08 0.999999987
[20,] 9.936044e-09 1.987209e-08 0.999999990
[21,] 6.038187e-09 1.207637e-08 0.999999994
[22,] 2.921917e-09 5.843835e-09 0.999999997
[23,] 1.920067e-09 3.840134e-09 0.999999998
[24,] 1.958639e-09 3.917277e-09 0.999999998
[25,] 4.239703e-09 8.479405e-09 0.999999996
[26,] 1.176390e-08 2.352779e-08 0.999999988
[27,] 2.408965e-08 4.817930e-08 0.999999976
[28,] 3.210499e-08 6.420998e-08 0.999999968
[29,] 2.200407e-08 4.400813e-08 0.999999978
[30,] 7.419905e-08 1.483981e-07 0.999999926
[31,] 4.126810e-07 8.253621e-07 0.999999587
[32,] 3.987330e-06 7.974659e-06 0.999996013
[33,] 2.435833e-05 4.871667e-05 0.999975642
[34,] 2.365800e-04 4.731599e-04 0.999763420
[35,] 1.391403e-03 2.782807e-03 0.998608597
[36,] 4.196980e-03 8.393959e-03 0.995803020
[37,] 1.180588e-02 2.361176e-02 0.988194118
[38,] 3.467536e-02 6.935073e-02 0.965324636
[39,] 3.690193e-02 7.380387e-02 0.963098067
[40,] 6.424164e-02 1.284833e-01 0.935758359
[41,] 1.308236e-01 2.616473e-01 0.869176359
[42,] 2.518426e-01 5.036853e-01 0.748157358
[43,] 5.145940e-01 9.708120e-01 0.485405981
[44,] 7.493829e-01 5.012342e-01 0.250617088
[45,] 8.823042e-01 2.353917e-01 0.117695832
[46,] 9.845362e-01 3.092758e-02 0.015463792
[47,] 9.959571e-01 8.085732e-03 0.004042866
[48,] 9.973963e-01 5.207466e-03 0.002603733
> postscript(file="/var/www/html/rcomp/tmp/1z6gg1258561768.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/2k1bw1258561768.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/3ak9g1258561768.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/4qo7s1258561768.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/56oob1258561768.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 = 79
Frequency = 1
1 2 3 4 5 6
7196.23413 7697.63358 7932.07276 10336.06414 10260.38303 8580.19596
7 8 9 10 11 12
7056.31643 -2218.94390 3476.80009 740.90507 4361.45237 4235.44758
13 14 15 16 17 18
3151.37512 4068.05052 4991.98984 6640.91579 6996.23844 7103.54024
19 20 21 22 23 24
4832.53495 -2421.09331 2774.12721 13260.37406 12492.35470 13441.98495
25 26 27 28 29 30
11453.74570 12309.59219 12715.05336 12629.60254 12176.04478 12317.82661
31 32 33 34 35 36
11222.87791 3544.83025 6523.43988 3337.74197 3751.15142 5252.85514
37 38 39 40 41 42
7618.19488 6117.84024 4338.49852 3181.45477 2471.49529 3070.59731
43 44 45 46 47 48
4757.47270 1428.07106 -1792.56012 -989.53824 -1640.27521 372.06060
49 50 51 52 53 54
-1974.38219 -1286.75773 1287.31327 438.97427 -323.27884 3656.28580
55 56 57 58 59 60
5642.02173 2147.84106 -1648.71314 -3502.22548 -2533.25234 -3572.71167
61 62 63 64 65 66
-2938.75382 -2714.85663 -21.49735 -1348.63094 706.18086 349.65294
67 68 69 70 71 72
786.05206 -2480.70517 -9333.09392 -12847.25738 -16431.43095 -19729.63661
73 74 75 76 77 78
-24506.41382 -26191.50218 -31243.43041 -31878.38057 -32287.06356 -35078.09886
79
-34297.27577
> postscript(file="/var/www/html/rcomp/tmp/6bdrr1258561768.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 = 79
Frequency = 1
lag(myerror, k = 1) myerror
0 7196.23413 NA
1 7697.63358 7196.23413
2 7932.07276 7697.63358
3 10336.06414 7932.07276
4 10260.38303 10336.06414
5 8580.19596 10260.38303
6 7056.31643 8580.19596
7 -2218.94390 7056.31643
8 3476.80009 -2218.94390
9 740.90507 3476.80009
10 4361.45237 740.90507
11 4235.44758 4361.45237
12 3151.37512 4235.44758
13 4068.05052 3151.37512
14 4991.98984 4068.05052
15 6640.91579 4991.98984
16 6996.23844 6640.91579
17 7103.54024 6996.23844
18 4832.53495 7103.54024
19 -2421.09331 4832.53495
20 2774.12721 -2421.09331
21 13260.37406 2774.12721
22 12492.35470 13260.37406
23 13441.98495 12492.35470
24 11453.74570 13441.98495
25 12309.59219 11453.74570
26 12715.05336 12309.59219
27 12629.60254 12715.05336
28 12176.04478 12629.60254
29 12317.82661 12176.04478
30 11222.87791 12317.82661
31 3544.83025 11222.87791
32 6523.43988 3544.83025
33 3337.74197 6523.43988
34 3751.15142 3337.74197
35 5252.85514 3751.15142
36 7618.19488 5252.85514
37 6117.84024 7618.19488
38 4338.49852 6117.84024
39 3181.45477 4338.49852
40 2471.49529 3181.45477
41 3070.59731 2471.49529
42 4757.47270 3070.59731
43 1428.07106 4757.47270
44 -1792.56012 1428.07106
45 -989.53824 -1792.56012
46 -1640.27521 -989.53824
47 372.06060 -1640.27521
48 -1974.38219 372.06060
49 -1286.75773 -1974.38219
50 1287.31327 -1286.75773
51 438.97427 1287.31327
52 -323.27884 438.97427
53 3656.28580 -323.27884
54 5642.02173 3656.28580
55 2147.84106 5642.02173
56 -1648.71314 2147.84106
57 -3502.22548 -1648.71314
58 -2533.25234 -3502.22548
59 -3572.71167 -2533.25234
60 -2938.75382 -3572.71167
61 -2714.85663 -2938.75382
62 -21.49735 -2714.85663
63 -1348.63094 -21.49735
64 706.18086 -1348.63094
65 349.65294 706.18086
66 786.05206 349.65294
67 -2480.70517 786.05206
68 -9333.09392 -2480.70517
69 -12847.25738 -9333.09392
70 -16431.43095 -12847.25738
71 -19729.63661 -16431.43095
72 -24506.41382 -19729.63661
73 -26191.50218 -24506.41382
74 -31243.43041 -26191.50218
75 -31878.38057 -31243.43041
76 -32287.06356 -31878.38057
77 -35078.09886 -32287.06356
78 -34297.27577 -35078.09886
79 NA -34297.27577
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 7697.63358 7196.23413
[2,] 7932.07276 7697.63358
[3,] 10336.06414 7932.07276
[4,] 10260.38303 10336.06414
[5,] 8580.19596 10260.38303
[6,] 7056.31643 8580.19596
[7,] -2218.94390 7056.31643
[8,] 3476.80009 -2218.94390
[9,] 740.90507 3476.80009
[10,] 4361.45237 740.90507
[11,] 4235.44758 4361.45237
[12,] 3151.37512 4235.44758
[13,] 4068.05052 3151.37512
[14,] 4991.98984 4068.05052
[15,] 6640.91579 4991.98984
[16,] 6996.23844 6640.91579
[17,] 7103.54024 6996.23844
[18,] 4832.53495 7103.54024
[19,] -2421.09331 4832.53495
[20,] 2774.12721 -2421.09331
[21,] 13260.37406 2774.12721
[22,] 12492.35470 13260.37406
[23,] 13441.98495 12492.35470
[24,] 11453.74570 13441.98495
[25,] 12309.59219 11453.74570
[26,] 12715.05336 12309.59219
[27,] 12629.60254 12715.05336
[28,] 12176.04478 12629.60254
[29,] 12317.82661 12176.04478
[30,] 11222.87791 12317.82661
[31,] 3544.83025 11222.87791
[32,] 6523.43988 3544.83025
[33,] 3337.74197 6523.43988
[34,] 3751.15142 3337.74197
[35,] 5252.85514 3751.15142
[36,] 7618.19488 5252.85514
[37,] 6117.84024 7618.19488
[38,] 4338.49852 6117.84024
[39,] 3181.45477 4338.49852
[40,] 2471.49529 3181.45477
[41,] 3070.59731 2471.49529
[42,] 4757.47270 3070.59731
[43,] 1428.07106 4757.47270
[44,] -1792.56012 1428.07106
[45,] -989.53824 -1792.56012
[46,] -1640.27521 -989.53824
[47,] 372.06060 -1640.27521
[48,] -1974.38219 372.06060
[49,] -1286.75773 -1974.38219
[50,] 1287.31327 -1286.75773
[51,] 438.97427 1287.31327
[52,] -323.27884 438.97427
[53,] 3656.28580 -323.27884
[54,] 5642.02173 3656.28580
[55,] 2147.84106 5642.02173
[56,] -1648.71314 2147.84106
[57,] -3502.22548 -1648.71314
[58,] -2533.25234 -3502.22548
[59,] -3572.71167 -2533.25234
[60,] -2938.75382 -3572.71167
[61,] -2714.85663 -2938.75382
[62,] -21.49735 -2714.85663
[63,] -1348.63094 -21.49735
[64,] 706.18086 -1348.63094
[65,] 349.65294 706.18086
[66,] 786.05206 349.65294
[67,] -2480.70517 786.05206
[68,] -9333.09392 -2480.70517
[69,] -12847.25738 -9333.09392
[70,] -16431.43095 -12847.25738
[71,] -19729.63661 -16431.43095
[72,] -24506.41382 -19729.63661
[73,] -26191.50218 -24506.41382
[74,] -31243.43041 -26191.50218
[75,] -31878.38057 -31243.43041
[76,] -32287.06356 -31878.38057
[77,] -35078.09886 -32287.06356
[78,] -34297.27577 -35078.09886
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 7697.63358 7196.23413
2 7932.07276 7697.63358
3 10336.06414 7932.07276
4 10260.38303 10336.06414
5 8580.19596 10260.38303
6 7056.31643 8580.19596
7 -2218.94390 7056.31643
8 3476.80009 -2218.94390
9 740.90507 3476.80009
10 4361.45237 740.90507
11 4235.44758 4361.45237
12 3151.37512 4235.44758
13 4068.05052 3151.37512
14 4991.98984 4068.05052
15 6640.91579 4991.98984
16 6996.23844 6640.91579
17 7103.54024 6996.23844
18 4832.53495 7103.54024
19 -2421.09331 4832.53495
20 2774.12721 -2421.09331
21 13260.37406 2774.12721
22 12492.35470 13260.37406
23 13441.98495 12492.35470
24 11453.74570 13441.98495
25 12309.59219 11453.74570
26 12715.05336 12309.59219
27 12629.60254 12715.05336
28 12176.04478 12629.60254
29 12317.82661 12176.04478
30 11222.87791 12317.82661
31 3544.83025 11222.87791
32 6523.43988 3544.83025
33 3337.74197 6523.43988
34 3751.15142 3337.74197
35 5252.85514 3751.15142
36 7618.19488 5252.85514
37 6117.84024 7618.19488
38 4338.49852 6117.84024
39 3181.45477 4338.49852
40 2471.49529 3181.45477
41 3070.59731 2471.49529
42 4757.47270 3070.59731
43 1428.07106 4757.47270
44 -1792.56012 1428.07106
45 -989.53824 -1792.56012
46 -1640.27521 -989.53824
47 372.06060 -1640.27521
48 -1974.38219 372.06060
49 -1286.75773 -1974.38219
50 1287.31327 -1286.75773
51 438.97427 1287.31327
52 -323.27884 438.97427
53 3656.28580 -323.27884
54 5642.02173 3656.28580
55 2147.84106 5642.02173
56 -1648.71314 2147.84106
57 -3502.22548 -1648.71314
58 -2533.25234 -3502.22548
59 -3572.71167 -2533.25234
60 -2938.75382 -3572.71167
61 -2714.85663 -2938.75382
62 -21.49735 -2714.85663
63 -1348.63094 -21.49735
64 706.18086 -1348.63094
65 349.65294 706.18086
66 786.05206 349.65294
67 -2480.70517 786.05206
68 -9333.09392 -2480.70517
69 -12847.25738 -9333.09392
70 -16431.43095 -12847.25738
71 -19729.63661 -16431.43095
72 -24506.41382 -19729.63661
73 -26191.50218 -24506.41382
74 -31243.43041 -26191.50218
75 -31878.38057 -31243.43041
76 -32287.06356 -31878.38057
77 -35078.09886 -32287.06356
78 -34297.27577 -35078.09886
> 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/7uuyc1258561768.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/8p1ak1258561768.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/96j0p1258561768.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/107rhj1258561768.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/11cu6o1258561768.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/12kzt61258561768.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/13pbb11258561768.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/147uu31258561768.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/15srap1258561768.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/167czh1258561768.tab")
+ }
> system("convert tmp/1z6gg1258561768.ps tmp/1z6gg1258561768.png")
> system("convert tmp/2k1bw1258561768.ps tmp/2k1bw1258561768.png")
> system("convert tmp/3ak9g1258561768.ps tmp/3ak9g1258561768.png")
> system("convert tmp/4qo7s1258561768.ps tmp/4qo7s1258561768.png")
> system("convert tmp/56oob1258561768.ps tmp/56oob1258561768.png")
> system("convert tmp/6bdrr1258561768.ps tmp/6bdrr1258561768.png")
> system("convert tmp/7uuyc1258561768.ps tmp/7uuyc1258561768.png")
> system("convert tmp/8p1ak1258561768.ps tmp/8p1ak1258561768.png")
> system("convert tmp/96j0p1258561768.ps tmp/96j0p1258561768.png")
> system("convert tmp/107rhj1258561768.ps tmp/107rhj1258561768.png")
>
>
> proc.time()
user system elapsed
2.723 1.598 3.782