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(209465
+ ,555332
+ ,213587
+ ,216234
+ ,204045
+ ,543599
+ ,209465
+ ,213587
+ ,200237
+ ,536662
+ ,204045
+ ,209465
+ ,203666
+ ,542722
+ ,200237
+ ,204045
+ ,241476
+ ,593530
+ ,203666
+ ,200237
+ ,260307
+ ,610763
+ ,241476
+ ,203666
+ ,243324
+ ,612613
+ ,260307
+ ,241476
+ ,244460
+ ,611324
+ ,243324
+ ,260307
+ ,233575
+ ,594167
+ ,244460
+ ,243324
+ ,237217
+ ,595454
+ ,233575
+ ,244460
+ ,235243
+ ,590865
+ ,237217
+ ,233575
+ ,230354
+ ,589379
+ ,235243
+ ,237217
+ ,227184
+ ,584428
+ ,230354
+ ,235243
+ ,221678
+ ,573100
+ ,227184
+ ,230354
+ ,217142
+ ,567456
+ ,221678
+ ,227184
+ ,219452
+ ,569028
+ ,217142
+ ,221678
+ ,256446
+ ,620735
+ ,219452
+ ,217142
+ ,265845
+ ,628884
+ ,256446
+ ,219452
+ ,248624
+ ,628232
+ ,265845
+ ,256446
+ ,241114
+ ,612117
+ ,248624
+ ,265845
+ ,229245
+ ,595404
+ ,241114
+ ,248624
+ ,231805
+ ,597141
+ ,229245
+ ,241114
+ ,219277
+ ,593408
+ ,231805
+ ,229245
+ ,219313
+ ,590072
+ ,219277
+ ,231805
+ ,212610
+ ,579799
+ ,219313
+ ,219277
+ ,214771
+ ,574205
+ ,212610
+ ,219313
+ ,211142
+ ,572775
+ ,214771
+ ,212610
+ ,211457
+ ,572942
+ ,211142
+ ,214771
+ ,240048
+ ,619567
+ ,211457
+ ,211142
+ ,240636
+ ,625809
+ ,240048
+ ,211457
+ ,230580
+ ,619916
+ ,240636
+ ,240048
+ ,208795
+ ,587625
+ ,230580
+ ,240636
+ ,197922
+ ,565742
+ ,208795
+ ,230580
+ ,194596
+ ,557274
+ ,197922
+ ,208795
+ ,194581
+ ,560576
+ ,194596
+ ,197922
+ ,185686
+ ,548854
+ ,194581
+ ,194596
+ ,178106
+ ,531673
+ ,185686
+ ,194581
+ ,172608
+ ,525919
+ ,178106
+ ,185686
+ ,167302
+ ,511038
+ ,172608
+ ,178106
+ ,168053
+ ,498662
+ ,167302
+ ,172608
+ ,202300
+ ,555362
+ ,168053
+ ,167302
+ ,202388
+ ,564591
+ ,202300
+ ,168053
+ ,182516
+ ,541657
+ ,202388
+ ,202300
+ ,173476
+ ,527070
+ ,182516
+ ,202388
+ ,166444
+ ,509846
+ ,173476
+ ,182516
+ ,171297
+ ,514258
+ ,166444
+ ,173476
+ ,169701
+ ,516922
+ ,171297
+ ,166444
+ ,164182
+ ,507561
+ ,169701
+ ,171297
+ ,161914
+ ,492622
+ ,164182
+ ,169701
+ ,159612
+ ,490243
+ ,161914
+ ,164182
+ ,151001
+ ,469357
+ ,159612
+ ,161914
+ ,158114
+ ,477580
+ ,151001
+ ,159612
+ ,186530
+ ,528379
+ ,158114
+ ,151001
+ ,187069
+ ,533590
+ ,186530
+ ,158114
+ ,174330
+ ,517945
+ ,187069
+ ,186530
+ ,169362
+ ,506174
+ ,174330
+ ,187069
+ ,166827
+ ,501866
+ ,169362
+ ,174330
+ ,178037
+ ,516141
+ ,166827
+ ,169362
+ ,186412
+ ,528222
+ ,178037
+ ,166827
+ ,189226
+ ,532638
+ ,186412
+ ,178037
+ ,191563
+ ,536322
+ ,189226
+ ,186412
+ ,188906
+ ,536535
+ ,191563
+ ,189226
+ ,186005
+ ,523597
+ ,188906
+ ,191563
+ ,195309
+ ,536214
+ ,186005
+ ,188906
+ ,223532
+ ,586570
+ ,195309
+ ,186005
+ ,226899
+ ,596594
+ ,223532
+ ,195309
+ ,214126
+ ,580523
+ ,226899
+ ,223532)
+ ,dim=c(4
+ ,67)
+ ,dimnames=list(c('Werkl'
+ ,'x'
+ ,'y-1'
+ ,'y-2')
+ ,1:67))
> y <- array(NA,dim=c(4,67),dimnames=list(c('Werkl','x','y-1','y-2'),1:67))
> 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 y-1 y-2
1 209465 555332 213587 216234
2 204045 543599 209465 213587
3 200237 536662 204045 209465
4 203666 542722 200237 204045
5 241476 593530 203666 200237
6 260307 610763 241476 203666
7 243324 612613 260307 241476
8 244460 611324 243324 260307
9 233575 594167 244460 243324
10 237217 595454 233575 244460
11 235243 590865 237217 233575
12 230354 589379 235243 237217
13 227184 584428 230354 235243
14 221678 573100 227184 230354
15 217142 567456 221678 227184
16 219452 569028 217142 221678
17 256446 620735 219452 217142
18 265845 628884 256446 219452
19 248624 628232 265845 256446
20 241114 612117 248624 265845
21 229245 595404 241114 248624
22 231805 597141 229245 241114
23 219277 593408 231805 229245
24 219313 590072 219277 231805
25 212610 579799 219313 219277
26 214771 574205 212610 219313
27 211142 572775 214771 212610
28 211457 572942 211142 214771
29 240048 619567 211457 211142
30 240636 625809 240048 211457
31 230580 619916 240636 240048
32 208795 587625 230580 240636
33 197922 565742 208795 230580
34 194596 557274 197922 208795
35 194581 560576 194596 197922
36 185686 548854 194581 194596
37 178106 531673 185686 194581
38 172608 525919 178106 185686
39 167302 511038 172608 178106
40 168053 498662 167302 172608
41 202300 555362 168053 167302
42 202388 564591 202300 168053
43 182516 541657 202388 202300
44 173476 527070 182516 202388
45 166444 509846 173476 182516
46 171297 514258 166444 173476
47 169701 516922 171297 166444
48 164182 507561 169701 171297
49 161914 492622 164182 169701
50 159612 490243 161914 164182
51 151001 469357 159612 161914
52 158114 477580 151001 159612
53 186530 528379 158114 151001
54 187069 533590 186530 158114
55 174330 517945 187069 186530
56 169362 506174 174330 187069
57 166827 501866 169362 174330
58 178037 516141 166827 169362
59 186412 528222 178037 166827
60 189226 532638 186412 178037
61 191563 536322 189226 186412
62 188906 536535 191563 189226
63 186005 523597 188906 191563
64 195309 536214 186005 188906
65 223532 586570 195309 186005
66 226899 596594 223532 195309
67 214126 580523 226899 223532
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) x `y-1` `y-2`
-1.479e+05 5.804e-01 2.779e-01 -1.403e-01
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-14721 -6133 1046 4897 16345
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.479e+05 1.822e+04 -8.115 2.24e-11 ***
x 5.804e-01 5.693e-02 10.195 5.86e-15 ***
`y-1` 2.779e-01 1.156e-01 2.405 0.0191 *
`y-2` -1.403e-01 7.618e-02 -1.841 0.0703 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7360 on 63 degrees of freedom
Multiple R-squared: 0.9384, Adjusted R-squared: 0.9354
F-statistic: 319.7 on 3 and 63 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.0124535393 2.490708e-02 9.875465e-01
[2,] 0.1334829423 2.669659e-01 8.665171e-01
[3,] 0.0621035277 1.242071e-01 9.378965e-01
[4,] 0.0412308039 8.246161e-02 9.587692e-01
[5,] 0.0214595476 4.291910e-02 9.785405e-01
[6,] 0.0107329599 2.146592e-02 9.892670e-01
[7,] 0.0052561518 1.051230e-02 9.947438e-01
[8,] 0.0025370042 5.074008e-03 9.974630e-01
[9,] 0.0012269531 2.453906e-03 9.987730e-01
[10,] 0.0006800061 1.360012e-03 9.993200e-01
[11,] 0.0010241766 2.048353e-03 9.989758e-01
[12,] 0.0013018520 2.603704e-03 9.986981e-01
[13,] 0.0068923202 1.378464e-02 9.931077e-01
[14,] 0.0062748032 1.254961e-02 9.937252e-01
[15,] 0.0077783859 1.555677e-02 9.922216e-01
[16,] 0.0143733247 2.874665e-02 9.856267e-01
[17,] 0.4819732831 9.639466e-01 5.180267e-01
[18,] 0.6872997280 6.254005e-01 3.127003e-01
[19,] 0.8675719989 2.648560e-01 1.324280e-01
[20,] 0.9016512203 1.966976e-01 9.834878e-02
[21,] 0.9407006712 1.185987e-01 5.929933e-02
[22,] 0.9564541172 8.709177e-02 4.354588e-02
[23,] 0.9663444349 6.731113e-02 3.365557e-02
[24,] 0.9885755170 2.284897e-02 1.142448e-02
[25,] 0.9940785131 1.184297e-02 5.921487e-03
[26,] 0.9969981518 6.003696e-03 3.001848e-03
[27,] 0.9966827211 6.634558e-03 3.317279e-03
[28,] 0.9964026166 7.194767e-03 3.597383e-03
[29,] 0.9975402829 4.919434e-03 2.459717e-03
[30,] 0.9990828380 1.834324e-03 9.171620e-04
[31,] 0.9989862545 2.027491e-03 1.013745e-03
[32,] 0.9994066248 1.186750e-03 5.933752e-04
[33,] 0.9992401517 1.519697e-03 7.598483e-04
[34,] 0.9987828685 2.434263e-03 1.217131e-03
[35,] 0.9978187734 4.362453e-03 2.181227e-03
[36,] 0.9983430228 3.313954e-03 1.656977e-03
[37,] 0.9991431111 1.713778e-03 8.568889e-04
[38,] 0.9994879067 1.024187e-03 5.120933e-04
[39,] 0.9995289556 9.420887e-04 4.710444e-04
[40,] 0.9994496847 1.100631e-03 5.503153e-04
[41,] 0.9998513002 2.973996e-04 1.486998e-04
[42,] 0.9999860267 2.794653e-05 1.397326e-05
[43,] 0.9999654555 6.908892e-05 3.454446e-05
[44,] 0.9999372138 1.255724e-04 6.278620e-05
[45,] 0.9998663268 2.673465e-04 1.336732e-04
[46,] 0.9998336654 3.326692e-04 1.663346e-04
[47,] 0.9995848682 8.302637e-04 4.151318e-04
[48,] 0.9988053041 2.389392e-03 1.194696e-03
[49,] 0.9982752280 3.449544e-03 1.724772e-03
[50,] 0.9965847736 6.830453e-03 3.415226e-03
[51,] 0.9957137948 8.572410e-03 4.286205e-03
[52,] 0.9937591555 1.248169e-02 6.240844e-03
[53,] 0.9818930559 3.621389e-02 1.810694e-02
[54,] 0.9433901483 1.132197e-01 5.660985e-02
> postscript(file="/var/www/html/rcomp/tmp/1uixq1259925714.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/2zbsv1259925714.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/34o1f1259925714.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/4asa31259925714.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/506ir1259925714.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 = 67
Frequency = 1
1 2 3 4 5 6
5990.1081 8153.9273 9300.0866 9509.9875 16345.0185 15147.0391
7 8 9 10 11 12
-2838.9667 6406.9524 2781.2071 8860.8652 7010.9541 4043.9530
13 14 15 16 17 18
4829.2476 6092.9002 5918.0903 7803.9840 13510.1268 8222.0866
19 20 21 22 23 24
-6042.8889 1904.6230 -593.2426 3203.7871 -9534.2534 -3721.1122
25 26 27 28 29 30
-6229.4577 1046.1588 -3293.8611 -1764.0228 -829.7036 -11766.4246
31 32 33 34 35 36
-14554.6859 -14721.4466 -8250.1925 -6695.8685 -9228.2480 -11782.5146
37 38 39 40 41 42
-6921.0221 -8220.7252 -4425.5013 4211.6009 4598.1923 -10082.8995
43 44 45 46 47 48
-11864.5373 -6903.2872 -4214.2590 -1235.7091 -6713.1220 -5674.8294
49 50 51 52 53 54
2037.3825 972.1794 4804.5217 9215.3734 4963.9370 -4421.0942
55 56 57 58 59 60
-4243.4622 1236.2678 795.1251 3727.8480 1620.1346 1116.1981
61 62 63 64 65 66
1707.9282 -1327.4334 4346.7767 6761.6920 2766.4751 -6222.8920
67
-6645.0727
> postscript(file="/var/www/html/rcomp/tmp/6br341259925714.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 = 67
Frequency = 1
lag(myerror, k = 1) myerror
0 5990.1081 NA
1 8153.9273 5990.1081
2 9300.0866 8153.9273
3 9509.9875 9300.0866
4 16345.0185 9509.9875
5 15147.0391 16345.0185
6 -2838.9667 15147.0391
7 6406.9524 -2838.9667
8 2781.2071 6406.9524
9 8860.8652 2781.2071
10 7010.9541 8860.8652
11 4043.9530 7010.9541
12 4829.2476 4043.9530
13 6092.9002 4829.2476
14 5918.0903 6092.9002
15 7803.9840 5918.0903
16 13510.1268 7803.9840
17 8222.0866 13510.1268
18 -6042.8889 8222.0866
19 1904.6230 -6042.8889
20 -593.2426 1904.6230
21 3203.7871 -593.2426
22 -9534.2534 3203.7871
23 -3721.1122 -9534.2534
24 -6229.4577 -3721.1122
25 1046.1588 -6229.4577
26 -3293.8611 1046.1588
27 -1764.0228 -3293.8611
28 -829.7036 -1764.0228
29 -11766.4246 -829.7036
30 -14554.6859 -11766.4246
31 -14721.4466 -14554.6859
32 -8250.1925 -14721.4466
33 -6695.8685 -8250.1925
34 -9228.2480 -6695.8685
35 -11782.5146 -9228.2480
36 -6921.0221 -11782.5146
37 -8220.7252 -6921.0221
38 -4425.5013 -8220.7252
39 4211.6009 -4425.5013
40 4598.1923 4211.6009
41 -10082.8995 4598.1923
42 -11864.5373 -10082.8995
43 -6903.2872 -11864.5373
44 -4214.2590 -6903.2872
45 -1235.7091 -4214.2590
46 -6713.1220 -1235.7091
47 -5674.8294 -6713.1220
48 2037.3825 -5674.8294
49 972.1794 2037.3825
50 4804.5217 972.1794
51 9215.3734 4804.5217
52 4963.9370 9215.3734
53 -4421.0942 4963.9370
54 -4243.4622 -4421.0942
55 1236.2678 -4243.4622
56 795.1251 1236.2678
57 3727.8480 795.1251
58 1620.1346 3727.8480
59 1116.1981 1620.1346
60 1707.9282 1116.1981
61 -1327.4334 1707.9282
62 4346.7767 -1327.4334
63 6761.6920 4346.7767
64 2766.4751 6761.6920
65 -6222.8920 2766.4751
66 -6645.0727 -6222.8920
67 NA -6645.0727
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 8153.9273 5990.1081
[2,] 9300.0866 8153.9273
[3,] 9509.9875 9300.0866
[4,] 16345.0185 9509.9875
[5,] 15147.0391 16345.0185
[6,] -2838.9667 15147.0391
[7,] 6406.9524 -2838.9667
[8,] 2781.2071 6406.9524
[9,] 8860.8652 2781.2071
[10,] 7010.9541 8860.8652
[11,] 4043.9530 7010.9541
[12,] 4829.2476 4043.9530
[13,] 6092.9002 4829.2476
[14,] 5918.0903 6092.9002
[15,] 7803.9840 5918.0903
[16,] 13510.1268 7803.9840
[17,] 8222.0866 13510.1268
[18,] -6042.8889 8222.0866
[19,] 1904.6230 -6042.8889
[20,] -593.2426 1904.6230
[21,] 3203.7871 -593.2426
[22,] -9534.2534 3203.7871
[23,] -3721.1122 -9534.2534
[24,] -6229.4577 -3721.1122
[25,] 1046.1588 -6229.4577
[26,] -3293.8611 1046.1588
[27,] -1764.0228 -3293.8611
[28,] -829.7036 -1764.0228
[29,] -11766.4246 -829.7036
[30,] -14554.6859 -11766.4246
[31,] -14721.4466 -14554.6859
[32,] -8250.1925 -14721.4466
[33,] -6695.8685 -8250.1925
[34,] -9228.2480 -6695.8685
[35,] -11782.5146 -9228.2480
[36,] -6921.0221 -11782.5146
[37,] -8220.7252 -6921.0221
[38,] -4425.5013 -8220.7252
[39,] 4211.6009 -4425.5013
[40,] 4598.1923 4211.6009
[41,] -10082.8995 4598.1923
[42,] -11864.5373 -10082.8995
[43,] -6903.2872 -11864.5373
[44,] -4214.2590 -6903.2872
[45,] -1235.7091 -4214.2590
[46,] -6713.1220 -1235.7091
[47,] -5674.8294 -6713.1220
[48,] 2037.3825 -5674.8294
[49,] 972.1794 2037.3825
[50,] 4804.5217 972.1794
[51,] 9215.3734 4804.5217
[52,] 4963.9370 9215.3734
[53,] -4421.0942 4963.9370
[54,] -4243.4622 -4421.0942
[55,] 1236.2678 -4243.4622
[56,] 795.1251 1236.2678
[57,] 3727.8480 795.1251
[58,] 1620.1346 3727.8480
[59,] 1116.1981 1620.1346
[60,] 1707.9282 1116.1981
[61,] -1327.4334 1707.9282
[62,] 4346.7767 -1327.4334
[63,] 6761.6920 4346.7767
[64,] 2766.4751 6761.6920
[65,] -6222.8920 2766.4751
[66,] -6645.0727 -6222.8920
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 8153.9273 5990.1081
2 9300.0866 8153.9273
3 9509.9875 9300.0866
4 16345.0185 9509.9875
5 15147.0391 16345.0185
6 -2838.9667 15147.0391
7 6406.9524 -2838.9667
8 2781.2071 6406.9524
9 8860.8652 2781.2071
10 7010.9541 8860.8652
11 4043.9530 7010.9541
12 4829.2476 4043.9530
13 6092.9002 4829.2476
14 5918.0903 6092.9002
15 7803.9840 5918.0903
16 13510.1268 7803.9840
17 8222.0866 13510.1268
18 -6042.8889 8222.0866
19 1904.6230 -6042.8889
20 -593.2426 1904.6230
21 3203.7871 -593.2426
22 -9534.2534 3203.7871
23 -3721.1122 -9534.2534
24 -6229.4577 -3721.1122
25 1046.1588 -6229.4577
26 -3293.8611 1046.1588
27 -1764.0228 -3293.8611
28 -829.7036 -1764.0228
29 -11766.4246 -829.7036
30 -14554.6859 -11766.4246
31 -14721.4466 -14554.6859
32 -8250.1925 -14721.4466
33 -6695.8685 -8250.1925
34 -9228.2480 -6695.8685
35 -11782.5146 -9228.2480
36 -6921.0221 -11782.5146
37 -8220.7252 -6921.0221
38 -4425.5013 -8220.7252
39 4211.6009 -4425.5013
40 4598.1923 4211.6009
41 -10082.8995 4598.1923
42 -11864.5373 -10082.8995
43 -6903.2872 -11864.5373
44 -4214.2590 -6903.2872
45 -1235.7091 -4214.2590
46 -6713.1220 -1235.7091
47 -5674.8294 -6713.1220
48 2037.3825 -5674.8294
49 972.1794 2037.3825
50 4804.5217 972.1794
51 9215.3734 4804.5217
52 4963.9370 9215.3734
53 -4421.0942 4963.9370
54 -4243.4622 -4421.0942
55 1236.2678 -4243.4622
56 795.1251 1236.2678
57 3727.8480 795.1251
58 1620.1346 3727.8480
59 1116.1981 1620.1346
60 1707.9282 1116.1981
61 -1327.4334 1707.9282
62 4346.7767 -1327.4334
63 6761.6920 4346.7767
64 2766.4751 6761.6920
65 -6222.8920 2766.4751
66 -6645.0727 -6222.8920
> 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/78hmf1259925714.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/8znqx1259925714.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/9gnyk1259925714.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/10djxv1259925714.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/11n9li1259925714.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/1282jw1259925714.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/132v3i1259925714.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/14gbky1259925714.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/15xf4m1259925714.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/16wh7o1259925714.tab")
+ }
>
> system("convert tmp/1uixq1259925714.ps tmp/1uixq1259925714.png")
> system("convert tmp/2zbsv1259925714.ps tmp/2zbsv1259925714.png")
> system("convert tmp/34o1f1259925714.ps tmp/34o1f1259925714.png")
> system("convert tmp/4asa31259925714.ps tmp/4asa31259925714.png")
> system("convert tmp/506ir1259925714.ps tmp/506ir1259925714.png")
> system("convert tmp/6br341259925714.ps tmp/6br341259925714.png")
> system("convert tmp/78hmf1259925714.ps tmp/78hmf1259925714.png")
> system("convert tmp/8znqx1259925714.ps tmp/8znqx1259925714.png")
> system("convert tmp/9gnyk1259925714.ps tmp/9gnyk1259925714.png")
> system("convert tmp/10djxv1259925714.ps tmp/10djxv1259925714.png")
>
>
> proc.time()
user system elapsed
2.561 1.553 6.412