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(216234
+ ,562325
+ ,213587
+ ,560854
+ ,209465
+ ,555332
+ ,204045
+ ,543599
+ ,200237
+ ,536662
+ ,203666
+ ,542722
+ ,241476
+ ,593530
+ ,260307
+ ,610763
+ ,243324
+ ,612613
+ ,244460
+ ,611324
+ ,233575
+ ,594167
+ ,237217
+ ,595454
+ ,235243
+ ,590865
+ ,230354
+ ,589379
+ ,227184
+ ,584428
+ ,221678
+ ,573100
+ ,217142
+ ,567456
+ ,219452
+ ,569028
+ ,256446
+ ,620735
+ ,265845
+ ,628884
+ ,248624
+ ,628232
+ ,241114
+ ,612117
+ ,229245
+ ,595404
+ ,231805
+ ,597141
+ ,219277
+ ,593408
+ ,219313
+ ,590072
+ ,212610
+ ,579799
+ ,214771
+ ,574205
+ ,211142
+ ,572775
+ ,211457
+ ,572942
+ ,240048
+ ,619567
+ ,240636
+ ,625809
+ ,230580
+ ,619916
+ ,208795
+ ,587625
+ ,197922
+ ,565742
+ ,194596
+ ,557274
+ ,194581
+ ,560576
+ ,185686
+ ,548854
+ ,178106
+ ,531673
+ ,172608
+ ,525919
+ ,167302
+ ,511038
+ ,168053
+ ,498662
+ ,202300
+ ,555362
+ ,202388
+ ,564591
+ ,182516
+ ,541657
+ ,173476
+ ,527070
+ ,166444
+ ,509846
+ ,171297
+ ,514258
+ ,169701
+ ,516922
+ ,164182
+ ,507561
+ ,161914
+ ,492622
+ ,159612
+ ,490243
+ ,151001
+ ,469357
+ ,158114
+ ,477580
+ ,186530
+ ,528379
+ ,187069
+ ,533590
+ ,174330
+ ,517945
+ ,169362
+ ,506174
+ ,166827
+ ,501866
+ ,178037
+ ,516141
+ ,186412
+ ,528222
+ ,189226
+ ,532638
+ ,191563
+ ,536322
+ ,188906
+ ,536535
+ ,186005
+ ,523597
+ ,195309
+ ,536214
+ ,223532
+ ,586570
+ ,226899
+ ,596594
+ ,214126
+ ,580523)
+ ,dim=c(2
+ ,69)
+ ,dimnames=list(c('Werkl'
+ ,'X')
+ ,1:69))
> y <- array(NA,dim=c(2,69),dimnames=list(c('Werkl','X'),1:69))
> 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 = 'Do not include Seasonal Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
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
Werkl X
1 216234 562325
2 213587 560854
3 209465 555332
4 204045 543599
5 200237 536662
6 203666 542722
7 241476 593530
8 260307 610763
9 243324 612613
10 244460 611324
11 233575 594167
12 237217 595454
13 235243 590865
14 230354 589379
15 227184 584428
16 221678 573100
17 217142 567456
18 219452 569028
19 256446 620735
20 265845 628884
21 248624 628232
22 241114 612117
23 229245 595404
24 231805 597141
25 219277 593408
26 219313 590072
27 212610 579799
28 214771 574205
29 211142 572775
30 211457 572942
31 240048 619567
32 240636 625809
33 230580 619916
34 208795 587625
35 197922 565742
36 194596 557274
37 194581 560576
38 185686 548854
39 178106 531673
40 172608 525919
41 167302 511038
42 168053 498662
43 202300 555362
44 202388 564591
45 182516 541657
46 173476 527070
47 166444 509846
48 171297 514258
49 169701 516922
50 164182 507561
51 161914 492622
52 159612 490243
53 151001 469357
54 158114 477580
55 186530 528379
56 187069 533590
57 174330 517945
58 169362 506174
59 166827 501866
60 178037 516141
61 186412 528222
62 189226 532638
63 191563 536322
64 188906 536535
65 186005 523597
66 195309 536214
67 223532 586570
68 226899 596594
69 214126 580523
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X
-1.764e+05 6.823e-01
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-16012.4 -5755.1 371.2 5823.1 19959.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.764e+05 1.276e+04 -13.82 <2e-16 ***
X 6.823e-01 2.278e-02 29.95 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7594 on 67 degrees of freedom
Multiple R-squared: 0.9305, Adjusted R-squared: 0.9295
F-statistic: 897 on 1 and 67 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.004919664 9.839328e-03 9.950803e-01
[2,] 0.000656724 1.313448e-03 9.993433e-01
[3,] 0.004437481 8.874962e-03 9.955625e-01
[4,] 0.016838385 3.367677e-02 9.831616e-01
[5,] 0.231858331 4.637167e-01 7.681417e-01
[6,] 0.233451714 4.669034e-01 7.665483e-01
[7,] 0.194739137 3.894783e-01 8.052609e-01
[8,] 0.144270666 2.885413e-01 8.557293e-01
[9,] 0.108804048 2.176081e-01 8.911960e-01
[10,] 0.089331326 1.786627e-01 9.106687e-01
[11,] 0.071798447 1.435969e-01 9.282016e-01
[12,] 0.054394794 1.087896e-01 9.456052e-01
[13,] 0.041884926 8.376985e-02 9.581151e-01
[14,] 0.033319209 6.663842e-02 9.666808e-01
[15,] 0.039257693 7.851539e-02 9.607423e-01
[16,] 0.133533198 2.670664e-01 8.664668e-01
[17,] 0.336648304 6.732966e-01 6.633517e-01
[18,] 0.415538203 8.310764e-01 5.844618e-01
[19,] 0.495624629 9.912493e-01 5.043754e-01
[20,] 0.561153155 8.776937e-01 4.388468e-01
[21,] 0.822426106 3.551478e-01 1.775739e-01
[22,] 0.894433227 2.111335e-01 1.055668e-01
[23,] 0.931167822 1.376644e-01 6.883218e-02
[24,] 0.932944856 1.341103e-01 6.705514e-02
[25,] 0.936618457 1.267631e-01 6.338154e-02
[26,] 0.937380886 1.252382e-01 6.261911e-02
[27,] 0.948591876 1.028162e-01 5.140812e-02
[28,] 0.960039236 7.992153e-02 3.996076e-02
[29,] 0.980909225 3.818155e-02 1.909077e-02
[30,] 0.994425387 1.114923e-02 5.574613e-03
[31,] 0.997220356 5.559288e-03 2.779644e-03
[32,] 0.997805041 4.389918e-03 2.194959e-03
[33,] 0.998786556 2.426888e-03 1.213444e-03
[34,] 0.999588837 8.223269e-04 4.111635e-04
[35,] 0.999667943 6.641132e-04 3.320566e-04
[36,] 0.999848794 3.024126e-04 1.512063e-04
[37,] 0.999810888 3.782246e-04 1.891123e-04
[38,] 0.999683239 6.335222e-04 3.167611e-04
[39,] 0.999409680 1.180640e-03 5.903200e-04
[40,] 0.999094953 1.810094e-03 9.050469e-04
[41,] 0.999664671 6.706588e-04 3.353294e-04
[42,] 0.999917424 1.651514e-04 8.257568e-05
[43,] 0.999931748 1.365036e-04 6.825180e-05
[44,] 0.999905859 1.882811e-04 9.414055e-05
[45,] 0.999973141 5.371791e-05 2.685896e-05
[46,] 0.999997019 5.961614e-06 2.980807e-06
[47,] 0.999992537 1.492637e-05 7.463184e-06
[48,] 0.999986028 2.794348e-05 1.397174e-05
[49,] 0.999965177 6.964511e-05 3.482256e-05
[50,] 0.999960223 7.955390e-05 3.977695e-05
[51,] 0.999885317 2.293670e-04 1.146835e-04
[52,] 0.999682440 6.351203e-04 3.175602e-04
[53,] 0.999722283 5.554346e-04 2.777173e-04
[54,] 0.999494975 1.010050e-03 5.050252e-04
[55,] 0.999411756 1.176487e-03 5.882436e-04
[56,] 0.998435694 3.128613e-03 1.564306e-03
[57,] 0.994942554 1.011489e-02 5.057446e-03
[58,] 0.984405659 3.118868e-02 1.559434e-02
[59,] 0.955038602 8.992280e-02 4.496140e-02
[60,] 0.944446682 1.111066e-01 5.555332e-02
> postscript(file="/var/www/html/rcomp/tmp/1awkz1259925096.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/21yqo1259925096.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/341381259925096.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/49j4s1259925096.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/57mkt1259925096.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 = 69
Frequency = 1
1 2 3 4 5 6
8933.6300 7290.2349 6935.6763 9520.6365 10445.4762 9739.9788
7 8 9 10 11 12
12885.6971 19959.3052 1714.1236 3729.5572 4550.0973 7314.0282
13 14 15 16 17 18
8470.9208 4595.7596 4803.6304 7026.2752 6340.9523 7578.4391
19 20 21 22 23 24
9294.8054 13134.0662 -3642.1001 -157.4753 -623.8587 751.0551
25 26 27 28 29 30
-9230.0672 -6918.0468 -6612.1866 -634.6225 -3287.9902 -3086.9277
31 32 33 34 35 36
-6306.3146 -9976.9834 -16012.4234 -15766.5558 -11709.6535 -9258.2732
37 38 39 40 41 42
-11526.0967 -12423.6414 -8281.7271 -9854.0013 -5007.2856 4187.3679
43 44 45 46 47 48
-249.7915 -6458.3719 -10683.4140 -9771.2829 -5052.0313 -3209.1638
49 50 51 52 53 54
-6622.7053 -5755.0665 2169.2202 1490.3175 7129.0062 8631.7797
55 56 57 58 59 60
2389.6384 -626.6201 -2691.6576 371.2285 775.4059 2246.1400
61 62 63 64 65 66
2378.7532 2179.8917 2003.4448 -798.8767 5127.2072 5823.1289
67 68 69
-309.7712 -3781.7485 -5590.1431
> postscript(file="/var/www/html/rcomp/tmp/6vueu1259925096.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 = 69
Frequency = 1
lag(myerror, k = 1) myerror
0 8933.6300 NA
1 7290.2349 8933.6300
2 6935.6763 7290.2349
3 9520.6365 6935.6763
4 10445.4762 9520.6365
5 9739.9788 10445.4762
6 12885.6971 9739.9788
7 19959.3052 12885.6971
8 1714.1236 19959.3052
9 3729.5572 1714.1236
10 4550.0973 3729.5572
11 7314.0282 4550.0973
12 8470.9208 7314.0282
13 4595.7596 8470.9208
14 4803.6304 4595.7596
15 7026.2752 4803.6304
16 6340.9523 7026.2752
17 7578.4391 6340.9523
18 9294.8054 7578.4391
19 13134.0662 9294.8054
20 -3642.1001 13134.0662
21 -157.4753 -3642.1001
22 -623.8587 -157.4753
23 751.0551 -623.8587
24 -9230.0672 751.0551
25 -6918.0468 -9230.0672
26 -6612.1866 -6918.0468
27 -634.6225 -6612.1866
28 -3287.9902 -634.6225
29 -3086.9277 -3287.9902
30 -6306.3146 -3086.9277
31 -9976.9834 -6306.3146
32 -16012.4234 -9976.9834
33 -15766.5558 -16012.4234
34 -11709.6535 -15766.5558
35 -9258.2732 -11709.6535
36 -11526.0967 -9258.2732
37 -12423.6414 -11526.0967
38 -8281.7271 -12423.6414
39 -9854.0013 -8281.7271
40 -5007.2856 -9854.0013
41 4187.3679 -5007.2856
42 -249.7915 4187.3679
43 -6458.3719 -249.7915
44 -10683.4140 -6458.3719
45 -9771.2829 -10683.4140
46 -5052.0313 -9771.2829
47 -3209.1638 -5052.0313
48 -6622.7053 -3209.1638
49 -5755.0665 -6622.7053
50 2169.2202 -5755.0665
51 1490.3175 2169.2202
52 7129.0062 1490.3175
53 8631.7797 7129.0062
54 2389.6384 8631.7797
55 -626.6201 2389.6384
56 -2691.6576 -626.6201
57 371.2285 -2691.6576
58 775.4059 371.2285
59 2246.1400 775.4059
60 2378.7532 2246.1400
61 2179.8917 2378.7532
62 2003.4448 2179.8917
63 -798.8767 2003.4448
64 5127.2072 -798.8767
65 5823.1289 5127.2072
66 -309.7712 5823.1289
67 -3781.7485 -309.7712
68 -5590.1431 -3781.7485
69 NA -5590.1431
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 7290.2349 8933.6300
[2,] 6935.6763 7290.2349
[3,] 9520.6365 6935.6763
[4,] 10445.4762 9520.6365
[5,] 9739.9788 10445.4762
[6,] 12885.6971 9739.9788
[7,] 19959.3052 12885.6971
[8,] 1714.1236 19959.3052
[9,] 3729.5572 1714.1236
[10,] 4550.0973 3729.5572
[11,] 7314.0282 4550.0973
[12,] 8470.9208 7314.0282
[13,] 4595.7596 8470.9208
[14,] 4803.6304 4595.7596
[15,] 7026.2752 4803.6304
[16,] 6340.9523 7026.2752
[17,] 7578.4391 6340.9523
[18,] 9294.8054 7578.4391
[19,] 13134.0662 9294.8054
[20,] -3642.1001 13134.0662
[21,] -157.4753 -3642.1001
[22,] -623.8587 -157.4753
[23,] 751.0551 -623.8587
[24,] -9230.0672 751.0551
[25,] -6918.0468 -9230.0672
[26,] -6612.1866 -6918.0468
[27,] -634.6225 -6612.1866
[28,] -3287.9902 -634.6225
[29,] -3086.9277 -3287.9902
[30,] -6306.3146 -3086.9277
[31,] -9976.9834 -6306.3146
[32,] -16012.4234 -9976.9834
[33,] -15766.5558 -16012.4234
[34,] -11709.6535 -15766.5558
[35,] -9258.2732 -11709.6535
[36,] -11526.0967 -9258.2732
[37,] -12423.6414 -11526.0967
[38,] -8281.7271 -12423.6414
[39,] -9854.0013 -8281.7271
[40,] -5007.2856 -9854.0013
[41,] 4187.3679 -5007.2856
[42,] -249.7915 4187.3679
[43,] -6458.3719 -249.7915
[44,] -10683.4140 -6458.3719
[45,] -9771.2829 -10683.4140
[46,] -5052.0313 -9771.2829
[47,] -3209.1638 -5052.0313
[48,] -6622.7053 -3209.1638
[49,] -5755.0665 -6622.7053
[50,] 2169.2202 -5755.0665
[51,] 1490.3175 2169.2202
[52,] 7129.0062 1490.3175
[53,] 8631.7797 7129.0062
[54,] 2389.6384 8631.7797
[55,] -626.6201 2389.6384
[56,] -2691.6576 -626.6201
[57,] 371.2285 -2691.6576
[58,] 775.4059 371.2285
[59,] 2246.1400 775.4059
[60,] 2378.7532 2246.1400
[61,] 2179.8917 2378.7532
[62,] 2003.4448 2179.8917
[63,] -798.8767 2003.4448
[64,] 5127.2072 -798.8767
[65,] 5823.1289 5127.2072
[66,] -309.7712 5823.1289
[67,] -3781.7485 -309.7712
[68,] -5590.1431 -3781.7485
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 7290.2349 8933.6300
2 6935.6763 7290.2349
3 9520.6365 6935.6763
4 10445.4762 9520.6365
5 9739.9788 10445.4762
6 12885.6971 9739.9788
7 19959.3052 12885.6971
8 1714.1236 19959.3052
9 3729.5572 1714.1236
10 4550.0973 3729.5572
11 7314.0282 4550.0973
12 8470.9208 7314.0282
13 4595.7596 8470.9208
14 4803.6304 4595.7596
15 7026.2752 4803.6304
16 6340.9523 7026.2752
17 7578.4391 6340.9523
18 9294.8054 7578.4391
19 13134.0662 9294.8054
20 -3642.1001 13134.0662
21 -157.4753 -3642.1001
22 -623.8587 -157.4753
23 751.0551 -623.8587
24 -9230.0672 751.0551
25 -6918.0468 -9230.0672
26 -6612.1866 -6918.0468
27 -634.6225 -6612.1866
28 -3287.9902 -634.6225
29 -3086.9277 -3287.9902
30 -6306.3146 -3086.9277
31 -9976.9834 -6306.3146
32 -16012.4234 -9976.9834
33 -15766.5558 -16012.4234
34 -11709.6535 -15766.5558
35 -9258.2732 -11709.6535
36 -11526.0967 -9258.2732
37 -12423.6414 -11526.0967
38 -8281.7271 -12423.6414
39 -9854.0013 -8281.7271
40 -5007.2856 -9854.0013
41 4187.3679 -5007.2856
42 -249.7915 4187.3679
43 -6458.3719 -249.7915
44 -10683.4140 -6458.3719
45 -9771.2829 -10683.4140
46 -5052.0313 -9771.2829
47 -3209.1638 -5052.0313
48 -6622.7053 -3209.1638
49 -5755.0665 -6622.7053
50 2169.2202 -5755.0665
51 1490.3175 2169.2202
52 7129.0062 1490.3175
53 8631.7797 7129.0062
54 2389.6384 8631.7797
55 -626.6201 2389.6384
56 -2691.6576 -626.6201
57 371.2285 -2691.6576
58 775.4059 371.2285
59 2246.1400 775.4059
60 2378.7532 2246.1400
61 2179.8917 2378.7532
62 2003.4448 2179.8917
63 -798.8767 2003.4448
64 5127.2072 -798.8767
65 5823.1289 5127.2072
66 -309.7712 5823.1289
67 -3781.7485 -309.7712
68 -5590.1431 -3781.7485
> 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/7td9e1259925096.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/8y56s1259925096.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/9bu1i1259925096.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/10pzto1259925096.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/118u1b1259925096.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/12w4re1259925096.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/133kpp1259925096.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/14ojdi1259925096.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/158r7t1259925096.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/16247c1259925096.tab")
+ }
>
> system("convert tmp/1awkz1259925096.ps tmp/1awkz1259925096.png")
> system("convert tmp/21yqo1259925096.ps tmp/21yqo1259925096.png")
> system("convert tmp/341381259925096.ps tmp/341381259925096.png")
> system("convert tmp/49j4s1259925096.ps tmp/49j4s1259925096.png")
> system("convert tmp/57mkt1259925096.ps tmp/57mkt1259925096.png")
> system("convert tmp/6vueu1259925096.ps tmp/6vueu1259925096.png")
> system("convert tmp/7td9e1259925096.ps tmp/7td9e1259925096.png")
> system("convert tmp/8y56s1259925096.ps tmp/8y56s1259925096.png")
> system("convert tmp/9bu1i1259925096.ps tmp/9bu1i1259925096.png")
> system("convert tmp/10pzto1259925096.ps tmp/10pzto1259925096.png")
>
>
> proc.time()
user system elapsed
2.523 1.563 3.221