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('Y'
+ ,'X')
+ ,1:69))
> y <- array(NA,dim=c(2,69),dimnames=list(c('Y','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 = '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 216234 562325 1 0 0 0 0 0 0 0 0 0 0
2 213587 560854 0 1 0 0 0 0 0 0 0 0 0
3 209465 555332 0 0 1 0 0 0 0 0 0 0 0
4 204045 543599 0 0 0 1 0 0 0 0 0 0 0
5 200237 536662 0 0 0 0 1 0 0 0 0 0 0
6 203666 542722 0 0 0 0 0 1 0 0 0 0 0
7 241476 593530 0 0 0 0 0 0 1 0 0 0 0
8 260307 610763 0 0 0 0 0 0 0 1 0 0 0
9 243324 612613 0 0 0 0 0 0 0 0 1 0 0
10 244460 611324 0 0 0 0 0 0 0 0 0 1 0
11 233575 594167 0 0 0 0 0 0 0 0 0 0 1
12 237217 595454 0 0 0 0 0 0 0 0 0 0 0
13 235243 590865 1 0 0 0 0 0 0 0 0 0 0
14 230354 589379 0 1 0 0 0 0 0 0 0 0 0
15 227184 584428 0 0 1 0 0 0 0 0 0 0 0
16 221678 573100 0 0 0 1 0 0 0 0 0 0 0
17 217142 567456 0 0 0 0 1 0 0 0 0 0 0
18 219452 569028 0 0 0 0 0 1 0 0 0 0 0
19 256446 620735 0 0 0 0 0 0 1 0 0 0 0
20 265845 628884 0 0 0 0 0 0 0 1 0 0 0
21 248624 628232 0 0 0 0 0 0 0 0 1 0 0
22 241114 612117 0 0 0 0 0 0 0 0 0 1 0
23 229245 595404 0 0 0 0 0 0 0 0 0 0 1
24 231805 597141 0 0 0 0 0 0 0 0 0 0 0
25 219277 593408 1 0 0 0 0 0 0 0 0 0 0
26 219313 590072 0 1 0 0 0 0 0 0 0 0 0
27 212610 579799 0 0 1 0 0 0 0 0 0 0 0
28 214771 574205 0 0 0 1 0 0 0 0 0 0 0
29 211142 572775 0 0 0 0 1 0 0 0 0 0 0
30 211457 572942 0 0 0 0 0 1 0 0 0 0 0
31 240048 619567 0 0 0 0 0 0 1 0 0 0 0
32 240636 625809 0 0 0 0 0 0 0 1 0 0 0
33 230580 619916 0 0 0 0 0 0 0 0 1 0 0
34 208795 587625 0 0 0 0 0 0 0 0 0 1 0
35 197922 565742 0 0 0 0 0 0 0 0 0 0 1
36 194596 557274 0 0 0 0 0 0 0 0 0 0 0
37 194581 560576 1 0 0 0 0 0 0 0 0 0 0
38 185686 548854 0 1 0 0 0 0 0 0 0 0 0
39 178106 531673 0 0 1 0 0 0 0 0 0 0 0
40 172608 525919 0 0 0 1 0 0 0 0 0 0 0
41 167302 511038 0 0 0 0 1 0 0 0 0 0 0
42 168053 498662 0 0 0 0 0 1 0 0 0 0 0
43 202300 555362 0 0 0 0 0 0 1 0 0 0 0
44 202388 564591 0 0 0 0 0 0 0 1 0 0 0
45 182516 541657 0 0 0 0 0 0 0 0 1 0 0
46 173476 527070 0 0 0 0 0 0 0 0 0 1 0
47 166444 509846 0 0 0 0 0 0 0 0 0 0 1
48 171297 514258 0 0 0 0 0 0 0 0 0 0 0
49 169701 516922 1 0 0 0 0 0 0 0 0 0 0
50 164182 507561 0 1 0 0 0 0 0 0 0 0 0
51 161914 492622 0 0 1 0 0 0 0 0 0 0 0
52 159612 490243 0 0 0 1 0 0 0 0 0 0 0
53 151001 469357 0 0 0 0 1 0 0 0 0 0 0
54 158114 477580 0 0 0 0 0 1 0 0 0 0 0
55 186530 528379 0 0 0 0 0 0 1 0 0 0 0
56 187069 533590 0 0 0 0 0 0 0 1 0 0 0
57 174330 517945 0 0 0 0 0 0 0 0 1 0 0
58 169362 506174 0 0 0 0 0 0 0 0 0 1 0
59 166827 501866 0 0 0 0 0 0 0 0 0 0 1
60 178037 516141 0 0 0 0 0 0 0 0 0 0 0
61 186412 528222 1 0 0 0 0 0 0 0 0 0 0
62 189226 532638 0 1 0 0 0 0 0 0 0 0 0
63 191563 536322 0 0 1 0 0 0 0 0 0 0 0
64 188906 536535 0 0 0 1 0 0 0 0 0 0 0
65 186005 523597 0 0 0 0 1 0 0 0 0 0 0
66 195309 536214 0 0 0 0 0 1 0 0 0 0 0
67 223532 586570 0 0 0 0 0 0 1 0 0 0 0
68 226899 596594 0 0 0 0 0 0 0 1 0 0 0
69 214126 580523 0 0 0 0 0 0 0 0 1 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
-1.882e+05 7.028e-01 -8.895e+02 -1.383e+03 7.934e+02 1.874e+03
M5 M6 M7 M8 M9 M10
4.422e+03 6.387e+03 2.806e+03 1.705e+03 -6.284e+03 -4.151e+03
M11
-1.926e+03
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-12686 -5803 924 4922 17560
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.882e+05 1.469e+04 -12.810 <2e-16 ***
X 7.028e-01 2.574e-02 27.310 <2e-16 ***
M1 -8.895e+02 4.514e+03 -0.197 0.844
M2 -1.383e+03 4.514e+03 -0.306 0.760
M3 7.934e+02 4.520e+03 0.176 0.861
M4 1.874e+03 4.531e+03 0.414 0.681
M5 4.422e+03 4.562e+03 0.969 0.337
M6 6.387e+03 4.553e+03 1.403 0.166
M7 2.806e+03 4.571e+03 0.614 0.542
M8 1.705e+03 4.615e+03 0.369 0.713
M9 -6.284e+03 4.568e+03 -1.376 0.174
M10 -4.151e+03 4.726e+03 -0.878 0.383
M11 -1.926e+03 4.715e+03 -0.409 0.684
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7454 on 56 degrees of freedom
Multiple R-squared: 0.944, Adjusted R-squared: 0.932
F-statistic: 78.72 on 12 and 56 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,] 2.207700e-03 4.415401e-03 9.977923e-01
[2,] 9.040803e-04 1.808161e-03 9.990959e-01
[3,] 1.135761e-04 2.271522e-04 9.998864e-01
[4,] 3.232893e-05 6.465785e-05 9.999677e-01
[5,] 1.464589e-03 2.929177e-03 9.985354e-01
[6,] 1.681908e-03 3.363816e-03 9.983181e-01
[7,] 2.728697e-03 5.457395e-03 9.972713e-01
[8,] 6.294510e-03 1.258902e-02 9.937055e-01
[9,] 2.121992e-02 4.243985e-02 9.787801e-01
[10,] 5.021319e-01 9.957362e-01 4.978681e-01
[11,] 6.440638e-01 7.118724e-01 3.559362e-01
[12,] 7.311408e-01 5.377184e-01 2.688592e-01
[13,] 7.628756e-01 4.742488e-01 2.371244e-01
[14,] 7.489287e-01 5.021425e-01 2.510713e-01
[15,] 7.370188e-01 5.259625e-01 2.629812e-01
[16,] 8.368797e-01 3.262406e-01 1.631203e-01
[17,] 9.760492e-01 4.790158e-02 2.395079e-02
[18,] 9.817206e-01 3.655878e-02 1.827939e-02
[19,] 9.946907e-01 1.061862e-02 5.309311e-03
[20,] 9.964463e-01 7.107377e-03 3.553689e-03
[21,] 9.971997e-01 5.600547e-03 2.800273e-03
[22,] 9.976816e-01 4.636734e-03 2.318367e-03
[23,] 9.984760e-01 3.047961e-03 1.523981e-03
[24,] 9.989032e-01 2.193582e-03 1.096791e-03
[25,] 9.994237e-01 1.152572e-03 5.762862e-04
[26,] 9.997757e-01 4.485419e-04 2.242710e-04
[27,] 9.995336e-01 9.328257e-04 4.664128e-04
[28,] 9.988419e-01 2.316280e-03 1.158140e-03
[29,] 9.980966e-01 3.806730e-03 1.903365e-03
[30,] 9.979074e-01 4.185284e-03 2.092642e-03
[31,] 9.986465e-01 2.706935e-03 1.353467e-03
[32,] 9.979045e-01 4.191077e-03 2.095538e-03
[33,] 9.967090e-01 6.582030e-03 3.291015e-03
[34,] 9.990377e-01 1.924500e-03 9.622502e-04
[35,] 9.999997e-01 6.951341e-07 3.475670e-07
[36,] 1.000000e+00 3.406849e-08 1.703424e-08
[37,] 9.999993e-01 1.402211e-06 7.011056e-07
[38,] 9.999989e-01 2.229824e-06 1.114912e-06
> postscript(file="/var/www/html/rcomp/tmp/1vz381258570572.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/2ptc71258570572.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/3ck7i1258570572.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/4oc271258570572.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/5yy1a1258570572.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
10125.3807 9006.0702 6588.3213 8333.9831 6853.5513 4058.0738
7 8 9 10 11 12
9739.4048 17559.8881 7265.6920 7174.8415 6123.4601 6934.6396
13 14 15 16 17 18
9075.4841 5724.7161 3857.6487 5232.6626 2115.4658 1355.3095
19 20 21 22 23 24
5588.7922 10361.8241 1588.1200 3271.4937 924.0539 336.9578
25 26 27 28 29 30
-8677.8240 -5803.3482 -7462.9304 -2450.9695 -7622.9109 -9390.5847
31 32 33 34 35 36
-9988.2971 -12685.9598 -10611.1083 -11833.6837 -9551.4692 -8852.1392
37 38 39 40 41 42
-10298.3618 -10460.9158 -8142.3173 -10676.9029 -8072.0175 -588.0381
43 44 45 46 47 48
-2610.8081 -7907.8374 -3671.9832 -4592.5406 -1743.8321 -1918.0104
49 50 51 52 53 54
-4496.8246 -2942.7708 3112.0727 1401.4207 4921.8273 4290.1188
55 56 57 58 59 60
583.7752 -1438.2651 4807.6284 5979.8891 4247.7872 3498.5522
61 62 63 64 65 66
4272.1456 4476.2484 2047.2051 -1840.1940 1804.0841 275.1208
67 68 69
-3312.8671 -5889.6499 621.6511
> postscript(file="/var/www/html/rcomp/tmp/6xuzy1258570572.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 10125.3807 NA
1 9006.0702 10125.3807
2 6588.3213 9006.0702
3 8333.9831 6588.3213
4 6853.5513 8333.9831
5 4058.0738 6853.5513
6 9739.4048 4058.0738
7 17559.8881 9739.4048
8 7265.6920 17559.8881
9 7174.8415 7265.6920
10 6123.4601 7174.8415
11 6934.6396 6123.4601
12 9075.4841 6934.6396
13 5724.7161 9075.4841
14 3857.6487 5724.7161
15 5232.6626 3857.6487
16 2115.4658 5232.6626
17 1355.3095 2115.4658
18 5588.7922 1355.3095
19 10361.8241 5588.7922
20 1588.1200 10361.8241
21 3271.4937 1588.1200
22 924.0539 3271.4937
23 336.9578 924.0539
24 -8677.8240 336.9578
25 -5803.3482 -8677.8240
26 -7462.9304 -5803.3482
27 -2450.9695 -7462.9304
28 -7622.9109 -2450.9695
29 -9390.5847 -7622.9109
30 -9988.2971 -9390.5847
31 -12685.9598 -9988.2971
32 -10611.1083 -12685.9598
33 -11833.6837 -10611.1083
34 -9551.4692 -11833.6837
35 -8852.1392 -9551.4692
36 -10298.3618 -8852.1392
37 -10460.9158 -10298.3618
38 -8142.3173 -10460.9158
39 -10676.9029 -8142.3173
40 -8072.0175 -10676.9029
41 -588.0381 -8072.0175
42 -2610.8081 -588.0381
43 -7907.8374 -2610.8081
44 -3671.9832 -7907.8374
45 -4592.5406 -3671.9832
46 -1743.8321 -4592.5406
47 -1918.0104 -1743.8321
48 -4496.8246 -1918.0104
49 -2942.7708 -4496.8246
50 3112.0727 -2942.7708
51 1401.4207 3112.0727
52 4921.8273 1401.4207
53 4290.1188 4921.8273
54 583.7752 4290.1188
55 -1438.2651 583.7752
56 4807.6284 -1438.2651
57 5979.8891 4807.6284
58 4247.7872 5979.8891
59 3498.5522 4247.7872
60 4272.1456 3498.5522
61 4476.2484 4272.1456
62 2047.2051 4476.2484
63 -1840.1940 2047.2051
64 1804.0841 -1840.1940
65 275.1208 1804.0841
66 -3312.8671 275.1208
67 -5889.6499 -3312.8671
68 621.6511 -5889.6499
69 NA 621.6511
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 9006.0702 10125.3807
[2,] 6588.3213 9006.0702
[3,] 8333.9831 6588.3213
[4,] 6853.5513 8333.9831
[5,] 4058.0738 6853.5513
[6,] 9739.4048 4058.0738
[7,] 17559.8881 9739.4048
[8,] 7265.6920 17559.8881
[9,] 7174.8415 7265.6920
[10,] 6123.4601 7174.8415
[11,] 6934.6396 6123.4601
[12,] 9075.4841 6934.6396
[13,] 5724.7161 9075.4841
[14,] 3857.6487 5724.7161
[15,] 5232.6626 3857.6487
[16,] 2115.4658 5232.6626
[17,] 1355.3095 2115.4658
[18,] 5588.7922 1355.3095
[19,] 10361.8241 5588.7922
[20,] 1588.1200 10361.8241
[21,] 3271.4937 1588.1200
[22,] 924.0539 3271.4937
[23,] 336.9578 924.0539
[24,] -8677.8240 336.9578
[25,] -5803.3482 -8677.8240
[26,] -7462.9304 -5803.3482
[27,] -2450.9695 -7462.9304
[28,] -7622.9109 -2450.9695
[29,] -9390.5847 -7622.9109
[30,] -9988.2971 -9390.5847
[31,] -12685.9598 -9988.2971
[32,] -10611.1083 -12685.9598
[33,] -11833.6837 -10611.1083
[34,] -9551.4692 -11833.6837
[35,] -8852.1392 -9551.4692
[36,] -10298.3618 -8852.1392
[37,] -10460.9158 -10298.3618
[38,] -8142.3173 -10460.9158
[39,] -10676.9029 -8142.3173
[40,] -8072.0175 -10676.9029
[41,] -588.0381 -8072.0175
[42,] -2610.8081 -588.0381
[43,] -7907.8374 -2610.8081
[44,] -3671.9832 -7907.8374
[45,] -4592.5406 -3671.9832
[46,] -1743.8321 -4592.5406
[47,] -1918.0104 -1743.8321
[48,] -4496.8246 -1918.0104
[49,] -2942.7708 -4496.8246
[50,] 3112.0727 -2942.7708
[51,] 1401.4207 3112.0727
[52,] 4921.8273 1401.4207
[53,] 4290.1188 4921.8273
[54,] 583.7752 4290.1188
[55,] -1438.2651 583.7752
[56,] 4807.6284 -1438.2651
[57,] 5979.8891 4807.6284
[58,] 4247.7872 5979.8891
[59,] 3498.5522 4247.7872
[60,] 4272.1456 3498.5522
[61,] 4476.2484 4272.1456
[62,] 2047.2051 4476.2484
[63,] -1840.1940 2047.2051
[64,] 1804.0841 -1840.1940
[65,] 275.1208 1804.0841
[66,] -3312.8671 275.1208
[67,] -5889.6499 -3312.8671
[68,] 621.6511 -5889.6499
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 9006.0702 10125.3807
2 6588.3213 9006.0702
3 8333.9831 6588.3213
4 6853.5513 8333.9831
5 4058.0738 6853.5513
6 9739.4048 4058.0738
7 17559.8881 9739.4048
8 7265.6920 17559.8881
9 7174.8415 7265.6920
10 6123.4601 7174.8415
11 6934.6396 6123.4601
12 9075.4841 6934.6396
13 5724.7161 9075.4841
14 3857.6487 5724.7161
15 5232.6626 3857.6487
16 2115.4658 5232.6626
17 1355.3095 2115.4658
18 5588.7922 1355.3095
19 10361.8241 5588.7922
20 1588.1200 10361.8241
21 3271.4937 1588.1200
22 924.0539 3271.4937
23 336.9578 924.0539
24 -8677.8240 336.9578
25 -5803.3482 -8677.8240
26 -7462.9304 -5803.3482
27 -2450.9695 -7462.9304
28 -7622.9109 -2450.9695
29 -9390.5847 -7622.9109
30 -9988.2971 -9390.5847
31 -12685.9598 -9988.2971
32 -10611.1083 -12685.9598
33 -11833.6837 -10611.1083
34 -9551.4692 -11833.6837
35 -8852.1392 -9551.4692
36 -10298.3618 -8852.1392
37 -10460.9158 -10298.3618
38 -8142.3173 -10460.9158
39 -10676.9029 -8142.3173
40 -8072.0175 -10676.9029
41 -588.0381 -8072.0175
42 -2610.8081 -588.0381
43 -7907.8374 -2610.8081
44 -3671.9832 -7907.8374
45 -4592.5406 -3671.9832
46 -1743.8321 -4592.5406
47 -1918.0104 -1743.8321
48 -4496.8246 -1918.0104
49 -2942.7708 -4496.8246
50 3112.0727 -2942.7708
51 1401.4207 3112.0727
52 4921.8273 1401.4207
53 4290.1188 4921.8273
54 583.7752 4290.1188
55 -1438.2651 583.7752
56 4807.6284 -1438.2651
57 5979.8891 4807.6284
58 4247.7872 5979.8891
59 3498.5522 4247.7872
60 4272.1456 3498.5522
61 4476.2484 4272.1456
62 2047.2051 4476.2484
63 -1840.1940 2047.2051
64 1804.0841 -1840.1940
65 275.1208 1804.0841
66 -3312.8671 275.1208
67 -5889.6499 -3312.8671
68 621.6511 -5889.6499
> 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/7blg01258570572.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/89nop1258570572.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/9f2th1258570572.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/1038gm1258570572.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/11xt0q1258570572.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/122qcr1258570573.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/13qoop1258570573.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/14rwhp1258570573.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/15y1jh1258570573.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/16t04o1258570573.tab")
+ }
>
> system("convert tmp/1vz381258570572.ps tmp/1vz381258570572.png")
> system("convert tmp/2ptc71258570572.ps tmp/2ptc71258570572.png")
> system("convert tmp/3ck7i1258570572.ps tmp/3ck7i1258570572.png")
> system("convert tmp/4oc271258570572.ps tmp/4oc271258570572.png")
> system("convert tmp/5yy1a1258570572.ps tmp/5yy1a1258570572.png")
> system("convert tmp/6xuzy1258570572.ps tmp/6xuzy1258570572.png")
> system("convert tmp/7blg01258570572.ps tmp/7blg01258570572.png")
> system("convert tmp/89nop1258570572.ps tmp/89nop1258570572.png")
> system("convert tmp/9f2th1258570572.ps tmp/9f2th1258570572.png")
> system("convert tmp/1038gm1258570572.ps tmp/1038gm1258570572.png")
>
>
> proc.time()
user system elapsed
2.488 1.562 2.944