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 = '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 274412 244752 1 0 0 0 0 0 0 0 0 0 0 1
2 272433 244576 0 1 0 0 0 0 0 0 0 0 0 2
3 268361 241572 0 0 1 0 0 0 0 0 0 0 0 3
4 268586 240541 0 0 0 1 0 0 0 0 0 0 0 4
5 264768 236089 0 0 0 0 1 0 0 0 0 0 0 5
6 269974 236997 0 0 0 0 0 1 0 0 0 0 0 6
7 304744 264579 0 0 0 0 0 0 1 0 0 0 0 7
8 309365 270349 0 0 0 0 0 0 0 1 0 0 0 8
9 308347 269645 0 0 0 0 0 0 0 0 1 0 0 9
10 298427 267037 0 0 0 0 0 0 0 0 0 1 0 10
11 289231 258113 0 0 0 0 0 0 0 0 0 0 1 11
12 291975 262813 0 0 0 0 0 0 0 0 0 0 0 12
13 294912 267413 1 0 0 0 0 0 0 0 0 0 0 13
14 293488 267366 0 1 0 0 0 0 0 0 0 0 0 14
15 290555 264777 0 0 1 0 0 0 0 0 0 0 0 15
16 284736 258863 0 0 0 1 0 0 0 0 0 0 0 16
17 281818 254844 0 0 0 0 1 0 0 0 0 0 0 17
18 287854 254868 0 0 0 0 0 1 0 0 0 0 0 18
19 316263 277267 0 0 0 0 0 0 1 0 0 0 0 19
20 325412 285351 0 0 0 0 0 0 0 1 0 0 0 20
21 326011 286602 0 0 0 0 0 0 0 0 1 0 0 21
22 328282 283042 0 0 0 0 0 0 0 0 0 1 0 22
23 317480 276687 0 0 0 0 0 0 0 0 0 0 1 23
24 317539 277915 0 0 0 0 0 0 0 0 0 0 0 24
25 313737 277128 1 0 0 0 0 0 0 0 0 0 0 25
26 312276 277103 0 1 0 0 0 0 0 0 0 0 0 26
27 309391 275037 0 0 1 0 0 0 0 0 0 0 0 27
28 302950 270150 0 0 0 1 0 0 0 0 0 0 0 28
29 300316 267140 0 0 0 0 1 0 0 0 0 0 0 29
30 304035 264993 0 0 0 0 0 1 0 0 0 0 0 30
31 333476 287259 0 0 0 0 0 0 1 0 0 0 0 31
32 337698 291186 0 0 0 0 0 0 0 1 0 0 0 32
33 335932 292300 0 0 0 0 0 0 0 0 1 0 0 33
34 323931 288186 0 0 0 0 0 0 0 0 0 1 0 34
35 313927 281477 0 0 0 0 0 0 0 0 0 0 1 35
36 314485 282656 0 0 0 0 0 0 0 0 0 0 0 36
37 313218 280190 1 0 0 0 0 0 0 0 0 0 0 37
38 309664 280408 0 1 0 0 0 0 0 0 0 0 0 38
39 302963 276836 0 0 1 0 0 0 0 0 0 0 0 39
40 298989 275216 0 0 0 1 0 0 0 0 0 0 0 40
41 298423 274352 0 0 0 0 1 0 0 0 0 0 0 41
42 301631 271311 0 0 0 0 0 1 0 0 0 0 0 42
43 329765 289802 0 0 0 0 0 0 1 0 0 0 0 43
44 335083 290726 0 0 0 0 0 0 0 1 0 0 0 44
45 327616 292300 0 0 0 0 0 0 0 0 1 0 0 45
46 309119 278506 0 0 0 0 0 0 0 0 0 1 0 46
47 295916 269826 0 0 0 0 0 0 0 0 0 0 1 47
48 291413 265861 0 0 0 0 0 0 0 0 0 0 0 48
49 291542 269034 1 0 0 0 0 0 0 0 0 0 0 49
50 284678 264176 0 1 0 0 0 0 0 0 0 0 0 50
51 276475 255198 0 0 1 0 0 0 0 0 0 0 0 51
52 272566 253353 0 0 0 1 0 0 0 0 0 0 0 52
53 264981 246057 0 0 0 0 1 0 0 0 0 0 0 53
54 263290 235372 0 0 0 0 0 1 0 0 0 0 0 54
55 296806 258556 0 0 0 0 0 0 1 0 0 0 0 55
56 303598 260993 0 0 0 0 0 0 0 1 0 0 0 56
57 286994 254663 0 0 0 0 0 0 0 0 1 0 0 57
58 276427 250643 0 0 0 0 0 0 0 0 0 1 0 58
59 266424 243422 0 0 0 0 0 0 0 0 0 0 1 59
60 267153 247105 0 0 0 0 0 0 0 0 0 0 0 60
61 268381 248541 1 0 0 0 0 0 0 0 0 0 0 61
62 262522 245039 0 1 0 0 0 0 0 0 0 0 0 62
63 255542 237080 0 0 1 0 0 0 0 0 0 0 0 63
64 253158 237085 0 0 0 1 0 0 0 0 0 0 0 64
65 243803 225554 0 0 0 0 1 0 0 0 0 0 0 65
66 250741 226839 0 0 0 0 0 1 0 0 0 0 0 66
67 280445 247934 0 0 0 0 0 0 1 0 0 0 0 67
68 285257 248333 0 0 0 0 0 0 0 1 0 0 0 68
69 270976 246969 0 0 0 0 0 0 0 0 1 0 0 69
70 261076 245098 0 0 0 0 0 0 0 0 0 1 0 70
71 255603 246263 0 0 0 0 0 0 0 0 0 0 1 71
72 260376 255765 0 0 0 0 0 0 0 0 0 0 0 72
73 263903 264319 1 0 0 0 0 0 0 0 0 0 0 73
74 264291 268347 0 1 0 0 0 0 0 0 0 0 0 74
75 263276 273046 0 0 1 0 0 0 0 0 0 0 0 75
76 262572 273963 0 0 0 1 0 0 0 0 0 0 0 76
77 256167 267430 0 0 0 0 1 0 0 0 0 0 0 77
78 264221 271993 0 0 0 0 0 1 0 0 0 0 0 78
79 293860 292710 0 0 0 0 0 0 1 0 0 0 0 79
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
41433.5106 0.9982 -2917.0703 -4883.4292 -5844.6724 -6705.0943
M5 M6 M7 M8 M9 M10
-5706.6546 462.0275 9147.7911 14952.8875 9314.8790 4907.5406
M11 t
1613.1280 -376.3602
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-19810 -5928 2705 6343 9433
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.143e+04 1.759e+04 2.355 0.02153 *
X 9.982e-01 6.361e-02 15.691 < 2e-16 ***
M1 -2.917e+03 4.749e+03 -0.614 0.54115
M2 -4.883e+03 4.748e+03 -1.029 0.30747
M3 -5.845e+03 4.756e+03 -1.229 0.22354
M4 -6.705e+03 4.765e+03 -1.407 0.16417
M5 -5.707e+03 4.808e+03 -1.187 0.23961
M6 4.620e+02 4.821e+03 0.096 0.92395
M7 9.148e+03 4.776e+03 1.915 0.05984 .
M8 1.495e+04 4.956e+03 3.017 0.00364 **
M9 9.315e+03 4.951e+03 1.882 0.06437 .
M10 4.908e+03 4.927e+03 0.996 0.32292
M11 1.613e+03 4.926e+03 0.328 0.74434
t -3.764e+02 4.269e+01 -8.817 1.04e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8525 on 65 degrees of freedom
Multiple R-squared: 0.8972, Adjusted R-squared: 0.8766
F-statistic: 43.64 on 13 and 65 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.0003031599 6.063198e-04 9.996968e-01
[2,] 0.0004135432 8.270864e-04 9.995865e-01
[3,] 0.0001004289 2.008578e-04 9.998996e-01
[4,] 0.0006683255 1.336651e-03 9.993317e-01
[5,] 0.0019124264 3.824853e-03 9.980876e-01
[6,] 0.8941329993 2.117340e-01 1.058670e-01
[7,] 0.9656856518 6.862870e-02 3.431435e-02
[8,] 0.9788523312 4.229534e-02 2.114767e-02
[9,] 0.9809428379 3.811432e-02 1.905716e-02
[10,] 0.9779959563 4.400809e-02 2.200404e-02
[11,] 0.9697225541 6.055489e-02 3.027745e-02
[12,] 0.9701197040 5.976059e-02 2.988030e-02
[13,] 0.9647407755 7.051845e-02 3.525922e-02
[14,] 0.9650766687 6.984666e-02 3.492333e-02
[15,] 0.9829062370 3.418753e-02 1.709376e-02
[16,] 0.9999079495 1.841011e-04 9.205054e-05
[17,] 0.9998929197 2.141606e-04 1.070803e-04
[18,] 0.9999906569 1.868620e-05 9.343102e-06
[19,] 0.9999972207 5.558667e-06 2.779333e-06
[20,] 0.9999972093 5.581389e-06 2.790694e-06
[21,] 0.9999953999 9.200116e-06 4.600058e-06
[22,] 0.9999931121 1.377583e-05 6.887913e-06
[23,] 0.9999900779 1.984410e-05 9.922051e-06
[24,] 0.9999915572 1.688554e-05 8.442770e-06
[25,] 0.9999930734 1.385320e-05 6.926602e-06
[26,] 0.9999896983 2.060348e-05 1.030174e-05
[27,] 0.9999734520 5.309591e-05 2.654795e-05
[28,] 0.9999504831 9.903389e-05 4.951695e-05
[29,] 0.9999936870 1.262600e-05 6.312998e-06
[30,] 0.9999993810 1.238050e-06 6.190249e-07
[31,] 0.9999998319 3.362732e-07 1.681366e-07
[32,] 0.9999999178 1.643763e-07 8.218815e-08
[33,] 0.9999998965 2.069809e-07 1.034905e-07
[34,] 0.9999996361 7.277066e-07 3.638533e-07
[35,] 0.9999988794 2.241143e-06 1.120571e-06
[36,] 0.9999956894 8.621242e-06 4.310621e-06
[37,] 0.9999895744 2.085113e-05 1.042556e-05
[38,] 0.9999999329 1.342456e-07 6.712279e-08
[39,] 0.9999999668 6.647834e-08 3.323917e-08
[40,] 0.9999999292 1.415553e-07 7.077767e-08
[41,] 0.9999995959 8.082063e-07 4.041032e-07
[42,] 0.9999965307 6.938510e-06 3.469255e-06
[43,] 0.9999707520 5.849595e-05 2.924798e-05
[44,] 0.9998788669 2.422662e-04 1.211331e-04
[45,] 0.9998034747 3.930506e-04 1.965253e-04
[46,] 0.9982000246 3.599951e-03 1.799975e-03
> postscript(file="/var/www/html/rcomp/tmp/155na1258903552.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/20rr31258903552.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/3ngq81258903552.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/4mw781258903552.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/5uthq1258903552.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
-8029.3072 -7489.9122 -7225.8414 -4734.9581 -4731.2371 -6223.8866
7 8 9 10 11 12
-7294.4907 -13861.5999 -8162.5278 -10695.6324 -7313.2955 -7271.1509
13 14 15 16 17 18
-5632.2483 -4666.6158 -3677.7806 -2356.8906 -1885.3721 -1665.6498
19 20 21 22 23 24
-3923.7997 -8272.6470 -2907.9743 7700.1677 6912.2361 7734.9861
25 26 27 28 29 30
8011.9672 8918.6402 9433.4386 9107.2201 8855.5969 8925.3208
31 32 33 34 35 36
7831.9259 2705.4224 5841.8429 2730.9645 3094.3810 4465.0407
37 38 39 40 41 42
8952.9295 7524.0501 5726.0748 4605.8733 4280.2025 4731.2799
43 44 45 46 47 48
6098.9324 5065.8981 2042.1657 2097.4592 1229.2456 2673.4313
49 50 51 52 53 54
2928.7057 3256.4776 5352.5458 4521.9299 3597.4130 6779.4120
55 56 57 58 59 60
8844.7078 7775.4601 3504.1705 1733.4653 2608.9388 1651.2102
61 62 63 64 65 66
4739.2855 4718.5548 7020.4997 5868.2911 7400.9743 7264.0191
67 68 69 70 71 72
7602.4675 6587.4663 -317.6771 -3566.4244 -6531.5060 -9253.5173
73 74 75 76 77 78
-10971.3324 -12261.1947 -16628.9369 -17011.4658 -17517.5775 -19810.4955
79
-19159.7432
> postscript(file="/var/www/html/rcomp/tmp/601uc1258903552.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 -8029.3072 NA
1 -7489.9122 -8029.3072
2 -7225.8414 -7489.9122
3 -4734.9581 -7225.8414
4 -4731.2371 -4734.9581
5 -6223.8866 -4731.2371
6 -7294.4907 -6223.8866
7 -13861.5999 -7294.4907
8 -8162.5278 -13861.5999
9 -10695.6324 -8162.5278
10 -7313.2955 -10695.6324
11 -7271.1509 -7313.2955
12 -5632.2483 -7271.1509
13 -4666.6158 -5632.2483
14 -3677.7806 -4666.6158
15 -2356.8906 -3677.7806
16 -1885.3721 -2356.8906
17 -1665.6498 -1885.3721
18 -3923.7997 -1665.6498
19 -8272.6470 -3923.7997
20 -2907.9743 -8272.6470
21 7700.1677 -2907.9743
22 6912.2361 7700.1677
23 7734.9861 6912.2361
24 8011.9672 7734.9861
25 8918.6402 8011.9672
26 9433.4386 8918.6402
27 9107.2201 9433.4386
28 8855.5969 9107.2201
29 8925.3208 8855.5969
30 7831.9259 8925.3208
31 2705.4224 7831.9259
32 5841.8429 2705.4224
33 2730.9645 5841.8429
34 3094.3810 2730.9645
35 4465.0407 3094.3810
36 8952.9295 4465.0407
37 7524.0501 8952.9295
38 5726.0748 7524.0501
39 4605.8733 5726.0748
40 4280.2025 4605.8733
41 4731.2799 4280.2025
42 6098.9324 4731.2799
43 5065.8981 6098.9324
44 2042.1657 5065.8981
45 2097.4592 2042.1657
46 1229.2456 2097.4592
47 2673.4313 1229.2456
48 2928.7057 2673.4313
49 3256.4776 2928.7057
50 5352.5458 3256.4776
51 4521.9299 5352.5458
52 3597.4130 4521.9299
53 6779.4120 3597.4130
54 8844.7078 6779.4120
55 7775.4601 8844.7078
56 3504.1705 7775.4601
57 1733.4653 3504.1705
58 2608.9388 1733.4653
59 1651.2102 2608.9388
60 4739.2855 1651.2102
61 4718.5548 4739.2855
62 7020.4997 4718.5548
63 5868.2911 7020.4997
64 7400.9743 5868.2911
65 7264.0191 7400.9743
66 7602.4675 7264.0191
67 6587.4663 7602.4675
68 -317.6771 6587.4663
69 -3566.4244 -317.6771
70 -6531.5060 -3566.4244
71 -9253.5173 -6531.5060
72 -10971.3324 -9253.5173
73 -12261.1947 -10971.3324
74 -16628.9369 -12261.1947
75 -17011.4658 -16628.9369
76 -17517.5775 -17011.4658
77 -19810.4955 -17517.5775
78 -19159.7432 -19810.4955
79 NA -19159.7432
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -7489.9122 -8029.3072
[2,] -7225.8414 -7489.9122
[3,] -4734.9581 -7225.8414
[4,] -4731.2371 -4734.9581
[5,] -6223.8866 -4731.2371
[6,] -7294.4907 -6223.8866
[7,] -13861.5999 -7294.4907
[8,] -8162.5278 -13861.5999
[9,] -10695.6324 -8162.5278
[10,] -7313.2955 -10695.6324
[11,] -7271.1509 -7313.2955
[12,] -5632.2483 -7271.1509
[13,] -4666.6158 -5632.2483
[14,] -3677.7806 -4666.6158
[15,] -2356.8906 -3677.7806
[16,] -1885.3721 -2356.8906
[17,] -1665.6498 -1885.3721
[18,] -3923.7997 -1665.6498
[19,] -8272.6470 -3923.7997
[20,] -2907.9743 -8272.6470
[21,] 7700.1677 -2907.9743
[22,] 6912.2361 7700.1677
[23,] 7734.9861 6912.2361
[24,] 8011.9672 7734.9861
[25,] 8918.6402 8011.9672
[26,] 9433.4386 8918.6402
[27,] 9107.2201 9433.4386
[28,] 8855.5969 9107.2201
[29,] 8925.3208 8855.5969
[30,] 7831.9259 8925.3208
[31,] 2705.4224 7831.9259
[32,] 5841.8429 2705.4224
[33,] 2730.9645 5841.8429
[34,] 3094.3810 2730.9645
[35,] 4465.0407 3094.3810
[36,] 8952.9295 4465.0407
[37,] 7524.0501 8952.9295
[38,] 5726.0748 7524.0501
[39,] 4605.8733 5726.0748
[40,] 4280.2025 4605.8733
[41,] 4731.2799 4280.2025
[42,] 6098.9324 4731.2799
[43,] 5065.8981 6098.9324
[44,] 2042.1657 5065.8981
[45,] 2097.4592 2042.1657
[46,] 1229.2456 2097.4592
[47,] 2673.4313 1229.2456
[48,] 2928.7057 2673.4313
[49,] 3256.4776 2928.7057
[50,] 5352.5458 3256.4776
[51,] 4521.9299 5352.5458
[52,] 3597.4130 4521.9299
[53,] 6779.4120 3597.4130
[54,] 8844.7078 6779.4120
[55,] 7775.4601 8844.7078
[56,] 3504.1705 7775.4601
[57,] 1733.4653 3504.1705
[58,] 2608.9388 1733.4653
[59,] 1651.2102 2608.9388
[60,] 4739.2855 1651.2102
[61,] 4718.5548 4739.2855
[62,] 7020.4997 4718.5548
[63,] 5868.2911 7020.4997
[64,] 7400.9743 5868.2911
[65,] 7264.0191 7400.9743
[66,] 7602.4675 7264.0191
[67,] 6587.4663 7602.4675
[68,] -317.6771 6587.4663
[69,] -3566.4244 -317.6771
[70,] -6531.5060 -3566.4244
[71,] -9253.5173 -6531.5060
[72,] -10971.3324 -9253.5173
[73,] -12261.1947 -10971.3324
[74,] -16628.9369 -12261.1947
[75,] -17011.4658 -16628.9369
[76,] -17517.5775 -17011.4658
[77,] -19810.4955 -17517.5775
[78,] -19159.7432 -19810.4955
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -7489.9122 -8029.3072
2 -7225.8414 -7489.9122
3 -4734.9581 -7225.8414
4 -4731.2371 -4734.9581
5 -6223.8866 -4731.2371
6 -7294.4907 -6223.8866
7 -13861.5999 -7294.4907
8 -8162.5278 -13861.5999
9 -10695.6324 -8162.5278
10 -7313.2955 -10695.6324
11 -7271.1509 -7313.2955
12 -5632.2483 -7271.1509
13 -4666.6158 -5632.2483
14 -3677.7806 -4666.6158
15 -2356.8906 -3677.7806
16 -1885.3721 -2356.8906
17 -1665.6498 -1885.3721
18 -3923.7997 -1665.6498
19 -8272.6470 -3923.7997
20 -2907.9743 -8272.6470
21 7700.1677 -2907.9743
22 6912.2361 7700.1677
23 7734.9861 6912.2361
24 8011.9672 7734.9861
25 8918.6402 8011.9672
26 9433.4386 8918.6402
27 9107.2201 9433.4386
28 8855.5969 9107.2201
29 8925.3208 8855.5969
30 7831.9259 8925.3208
31 2705.4224 7831.9259
32 5841.8429 2705.4224
33 2730.9645 5841.8429
34 3094.3810 2730.9645
35 4465.0407 3094.3810
36 8952.9295 4465.0407
37 7524.0501 8952.9295
38 5726.0748 7524.0501
39 4605.8733 5726.0748
40 4280.2025 4605.8733
41 4731.2799 4280.2025
42 6098.9324 4731.2799
43 5065.8981 6098.9324
44 2042.1657 5065.8981
45 2097.4592 2042.1657
46 1229.2456 2097.4592
47 2673.4313 1229.2456
48 2928.7057 2673.4313
49 3256.4776 2928.7057
50 5352.5458 3256.4776
51 4521.9299 5352.5458
52 3597.4130 4521.9299
53 6779.4120 3597.4130
54 8844.7078 6779.4120
55 7775.4601 8844.7078
56 3504.1705 7775.4601
57 1733.4653 3504.1705
58 2608.9388 1733.4653
59 1651.2102 2608.9388
60 4739.2855 1651.2102
61 4718.5548 4739.2855
62 7020.4997 4718.5548
63 5868.2911 7020.4997
64 7400.9743 5868.2911
65 7264.0191 7400.9743
66 7602.4675 7264.0191
67 6587.4663 7602.4675
68 -317.6771 6587.4663
69 -3566.4244 -317.6771
70 -6531.5060 -3566.4244
71 -9253.5173 -6531.5060
72 -10971.3324 -9253.5173
73 -12261.1947 -10971.3324
74 -16628.9369 -12261.1947
75 -17011.4658 -16628.9369
76 -17517.5775 -17011.4658
77 -19810.4955 -17517.5775
78 -19159.7432 -19810.4955
> 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/7uedv1258903552.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/8u8y41258903552.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/9fuga1258903552.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/109djl1258903552.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/11jfxr1258903552.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/12z3wo1258903552.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/13lhmg1258903552.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/14imff1258903552.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/15oa8t1258903552.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/164rqd1258903552.tab")
+ }
> system("convert tmp/155na1258903552.ps tmp/155na1258903552.png")
> system("convert tmp/20rr31258903552.ps tmp/20rr31258903552.png")
> system("convert tmp/3ngq81258903552.ps tmp/3ngq81258903552.png")
> system("convert tmp/4mw781258903552.ps tmp/4mw781258903552.png")
> system("convert tmp/5uthq1258903552.ps tmp/5uthq1258903552.png")
> system("convert tmp/601uc1258903552.ps tmp/601uc1258903552.png")
> system("convert tmp/7uedv1258903552.ps tmp/7uedv1258903552.png")
> system("convert tmp/8u8y41258903552.ps tmp/8u8y41258903552.png")
> system("convert tmp/9fuga1258903552.ps tmp/9fuga1258903552.png")
> system("convert tmp/109djl1258903552.ps tmp/109djl1258903552.png")
>
>
> proc.time()
user system elapsed
2.634 1.595 3.525