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(344744
+ ,492865
+ ,338653
+ ,480961
+ ,327532
+ ,461935
+ ,326225
+ ,456608
+ ,318672
+ ,441977
+ ,317756
+ ,439148
+ ,337302
+ ,488180
+ ,349420
+ ,520564
+ ,336923
+ ,501492
+ ,330758
+ ,485025
+ ,321002
+ ,464196
+ ,320820
+ ,460170
+ ,327032
+ ,467037
+ ,324047
+ ,460070
+ ,316735
+ ,447988
+ ,315710
+ ,442867
+ ,313427
+ ,436087
+ ,310527
+ ,431328
+ ,330962
+ ,484015
+ ,339015
+ ,509673
+ ,341332
+ ,512927
+ ,339092
+ ,502831
+ ,323308
+ ,470984
+ ,325849
+ ,471067
+ ,330675
+ ,476049
+ ,332225
+ ,474605
+ ,331735
+ ,470439
+ ,328047
+ ,461251
+ ,326165
+ ,454724
+ ,327081
+ ,455626
+ ,346764
+ ,516847
+ ,344190
+ ,525192
+ ,343333
+ ,522975
+ ,345777
+ ,518585
+ ,344094
+ ,509239
+ ,348609
+ ,512238
+ ,354846
+ ,519164
+ ,356427
+ ,517009
+ ,353467
+ ,509933
+ ,355996
+ ,509127
+ ,352487
+ ,500857
+ ,355178
+ ,506971
+ ,374556
+ ,569323
+ ,375021
+ ,579714
+ ,375787
+ ,577992
+ ,372720
+ ,565464
+ ,364431
+ ,547344
+ ,370490
+ ,554788
+ ,376974
+ ,562325
+ ,377632
+ ,560854
+ ,378205
+ ,555332
+ ,370861
+ ,543599
+ ,369167
+ ,536662
+ ,371551
+ ,542722
+ ,382842
+ ,593530
+ ,381903
+ ,610763
+ ,384502
+ ,612613
+ ,392058
+ ,611324
+ ,384359
+ ,594167
+ ,388884
+ ,595454
+ ,386586
+ ,590865
+ ,387495
+ ,589379
+ ,385705
+ ,584428
+ ,378670
+ ,573100
+ ,377367
+ ,567456
+ ,376911
+ ,569028
+ ,389827
+ ,620735
+ ,387820
+ ,628884
+ ,387267
+ ,628232
+ ,380575
+ ,612117
+ ,372402
+ ,595404
+ ,376740
+ ,597141)
+ ,dim=c(2
+ ,72)
+ ,dimnames=list(c('Y'
+ ,'X')
+ ,1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('Y','X'),1:72))
> 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'
> 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
1 344744 492865
2 338653 480961
3 327532 461935
4 326225 456608
5 318672 441977
6 317756 439148
7 337302 488180
8 349420 520564
9 336923 501492
10 330758 485025
11 321002 464196
12 320820 460170
13 327032 467037
14 324047 460070
15 316735 447988
16 315710 442867
17 313427 436087
18 310527 431328
19 330962 484015
20 339015 509673
21 341332 512927
22 339092 502831
23 323308 470984
24 325849 471067
25 330675 476049
26 332225 474605
27 331735 470439
28 328047 461251
29 326165 454724
30 327081 455626
31 346764 516847
32 344190 525192
33 343333 522975
34 345777 518585
35 344094 509239
36 348609 512238
37 354846 519164
38 356427 517009
39 353467 509933
40 355996 509127
41 352487 500857
42 355178 506971
43 374556 569323
44 375021 579714
45 375787 577992
46 372720 565464
47 364431 547344
48 370490 554788
49 376974 562325
50 377632 560854
51 378205 555332
52 370861 543599
53 369167 536662
54 371551 542722
55 382842 593530
56 381903 610763
57 384502 612613
58 392058 611324
59 384359 594167
60 388884 595454
61 386586 590865
62 387495 589379
63 385705 584428
64 378670 573100
65 377367 567456
66 376911 569028
67 389827 620735
68 387820 628884
69 387267 628232
70 380575 612117
71 372402 595404
72 376740 597141
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X
1.298e+05 4.260e-01
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-11055.7 -4937.8 398.3 4695.1 11817.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.298e+05 6.687e+03 19.42 <2e-16 ***
X 4.260e-01 1.268e-02 33.60 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6012 on 70 degrees of freedom
Multiple R-squared: 0.9416, Adjusted R-squared: 0.9408
F-statistic: 1129 on 1 and 70 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.0020802368 0.0041604736 0.99791976
[2,] 0.0003338901 0.0006677802 0.99966611
[3,] 0.0063354729 0.0126709457 0.99366453
[4,] 0.0116436823 0.0232873645 0.98835632
[5,] 0.0431254967 0.0862509935 0.95687450
[6,] 0.0536514147 0.1073028293 0.94634859
[7,] 0.0839274353 0.1678548705 0.91607256
[8,] 0.0759808605 0.1519617210 0.92401914
[9,] 0.0451890954 0.0903781907 0.95481090
[10,] 0.0259947167 0.0519894335 0.97400528
[11,] 0.0184957161 0.0369914322 0.98150428
[12,] 0.0109071677 0.0218143354 0.98909283
[13,] 0.0059501724 0.0119003447 0.99404983
[14,] 0.0034239390 0.0068478779 0.99657606
[15,] 0.0026563827 0.0053127655 0.99734362
[16,] 0.0038482381 0.0076964762 0.99615176
[17,] 0.0034529928 0.0069059855 0.99654701
[18,] 0.0022272911 0.0044545821 0.99777271
[19,] 0.0030843492 0.0061686983 0.99691565
[20,] 0.0024626305 0.0049252610 0.99753737
[21,] 0.0015900459 0.0031800917 0.99840995
[22,] 0.0011706473 0.0023412945 0.99882935
[23,] 0.0010092765 0.0020185529 0.99899072
[24,] 0.0008488303 0.0016976605 0.99915117
[25,] 0.0007903946 0.0015807892 0.99920961
[26,] 0.0008094807 0.0016189613 0.99919052
[27,] 0.0006342885 0.0012685769 0.99936571
[28,] 0.0019106009 0.0038212018 0.99808940
[29,] 0.0070653716 0.0141307432 0.99293463
[30,] 0.0129649244 0.0259298489 0.98703508
[31,] 0.0278951272 0.0557902545 0.97210487
[32,] 0.0586719259 0.1173438519 0.94132807
[33,] 0.1241073198 0.2482146396 0.87589268
[34,] 0.2392668549 0.4785337098 0.76073315
[35,] 0.3697660877 0.7395321753 0.63023391
[36,] 0.5349611222 0.9300777555 0.46503888
[37,] 0.6942391983 0.6115216034 0.30576080
[38,] 0.8164421712 0.3671156576 0.18355783
[39,] 0.7784162194 0.4431675612 0.22158378
[40,] 0.7491909950 0.5016180100 0.25080901
[41,] 0.7039118517 0.5921762966 0.29608815
[42,] 0.6677311558 0.6645376885 0.33226884
[43,] 0.7545692565 0.4908614870 0.24543074
[44,] 0.7533021278 0.4933957443 0.24669787
[45,] 0.7309031928 0.5381936144 0.26909681
[46,] 0.7161103678 0.5677792645 0.28388963
[47,] 0.7450496061 0.5099007879 0.25495039
[48,] 0.7242121853 0.5515756293 0.27578781
[49,] 0.7219478334 0.5561043332 0.27805217
[50,] 0.7048335456 0.5903329087 0.29516645
[51,] 0.6389487610 0.7221024779 0.36105124
[52,] 0.6680854398 0.6638291204 0.33191456
[53,] 0.6302633410 0.7394733180 0.36973666
[54,] 0.6522056925 0.6955886151 0.34779431
[55,] 0.5695423545 0.8609152911 0.43045765
[56,] 0.5976884056 0.8046231888 0.40231159
[57,] 0.5894580464 0.8210839073 0.41054195
[58,] 0.6650600690 0.6698798620 0.33493993
[59,] 0.7534198787 0.4931602427 0.24658012
[60,] 0.6831580357 0.6336839287 0.31684196
[61,] 0.6466184944 0.7067630112 0.35338151
[62,] 0.8992907900 0.2014184201 0.10070921
[63,] 0.9721258150 0.0557483699 0.02787418
> postscript(file="/var/www/html/rcomp/tmp/1zr841258484720.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/2lidu1258484720.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/391qq1258484720.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/48xiu1258484720.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/5pxk41258484720.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 = 72
Frequency = 1
1 2 3 4 5 6
4965.9573 3945.8338 929.5462 1891.7464 571.2728 860.3727
7 8 9 10 11 12
-480.3222 -2157.2873 -6529.9798 -5680.3525 -6563.5966 -5030.5975
13 14 15 16 17 18
-1743.8083 -1760.9994 -3926.2983 -2769.8502 -2164.6998 -3037.4567
19 20 21 22 23 24
-5046.1118 -7922.9295 -6992.0713 -4931.3682 -7149.1548 -4643.5113
25 26 27 28 29 30
-1939.7481 225.3684 1510.0047 1735.9171 2634.2944 3166.0596
31 32 33 34 35 36
-3229.9164 -9358.7269 -9271.3272 -4957.2711 -2659.0538 578.4296
37 38 39 40 41 42
3865.0859 6364.0747 6418.3155 9290.6561 9304.5181 9391.0709
43 44 45 46 47 48
2208.3105 -1753.0569 -253.5178 2016.1708 1445.9445 4333.9427
49 50 51 52 53 54
7607.3248 8891.9427 11817.2091 9471.2429 10732.2723 10534.8281
55 56 57 58 59 60
182.5911 -8097.3375 -6286.4022 1818.6872 1428.2413 5405.0039
61 62 63 64 65 66
5061.8302 6603.8378 6922.8692 4713.3807 5814.6168 4688.9749
67 68 69 70 71 72
-4421.2190 -9899.5372 -10174.7977 -10002.1157 -11055.6971 -7457.6259
> postscript(file="/var/www/html/rcomp/tmp/6fbjy1258484720.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 4965.9573 NA
1 3945.8338 4965.9573
2 929.5462 3945.8338
3 1891.7464 929.5462
4 571.2728 1891.7464
5 860.3727 571.2728
6 -480.3222 860.3727
7 -2157.2873 -480.3222
8 -6529.9798 -2157.2873
9 -5680.3525 -6529.9798
10 -6563.5966 -5680.3525
11 -5030.5975 -6563.5966
12 -1743.8083 -5030.5975
13 -1760.9994 -1743.8083
14 -3926.2983 -1760.9994
15 -2769.8502 -3926.2983
16 -2164.6998 -2769.8502
17 -3037.4567 -2164.6998
18 -5046.1118 -3037.4567
19 -7922.9295 -5046.1118
20 -6992.0713 -7922.9295
21 -4931.3682 -6992.0713
22 -7149.1548 -4931.3682
23 -4643.5113 -7149.1548
24 -1939.7481 -4643.5113
25 225.3684 -1939.7481
26 1510.0047 225.3684
27 1735.9171 1510.0047
28 2634.2944 1735.9171
29 3166.0596 2634.2944
30 -3229.9164 3166.0596
31 -9358.7269 -3229.9164
32 -9271.3272 -9358.7269
33 -4957.2711 -9271.3272
34 -2659.0538 -4957.2711
35 578.4296 -2659.0538
36 3865.0859 578.4296
37 6364.0747 3865.0859
38 6418.3155 6364.0747
39 9290.6561 6418.3155
40 9304.5181 9290.6561
41 9391.0709 9304.5181
42 2208.3105 9391.0709
43 -1753.0569 2208.3105
44 -253.5178 -1753.0569
45 2016.1708 -253.5178
46 1445.9445 2016.1708
47 4333.9427 1445.9445
48 7607.3248 4333.9427
49 8891.9427 7607.3248
50 11817.2091 8891.9427
51 9471.2429 11817.2091
52 10732.2723 9471.2429
53 10534.8281 10732.2723
54 182.5911 10534.8281
55 -8097.3375 182.5911
56 -6286.4022 -8097.3375
57 1818.6872 -6286.4022
58 1428.2413 1818.6872
59 5405.0039 1428.2413
60 5061.8302 5405.0039
61 6603.8378 5061.8302
62 6922.8692 6603.8378
63 4713.3807 6922.8692
64 5814.6168 4713.3807
65 4688.9749 5814.6168
66 -4421.2190 4688.9749
67 -9899.5372 -4421.2190
68 -10174.7977 -9899.5372
69 -10002.1157 -10174.7977
70 -11055.6971 -10002.1157
71 -7457.6259 -11055.6971
72 NA -7457.6259
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 3945.8338 4965.9573
[2,] 929.5462 3945.8338
[3,] 1891.7464 929.5462
[4,] 571.2728 1891.7464
[5,] 860.3727 571.2728
[6,] -480.3222 860.3727
[7,] -2157.2873 -480.3222
[8,] -6529.9798 -2157.2873
[9,] -5680.3525 -6529.9798
[10,] -6563.5966 -5680.3525
[11,] -5030.5975 -6563.5966
[12,] -1743.8083 -5030.5975
[13,] -1760.9994 -1743.8083
[14,] -3926.2983 -1760.9994
[15,] -2769.8502 -3926.2983
[16,] -2164.6998 -2769.8502
[17,] -3037.4567 -2164.6998
[18,] -5046.1118 -3037.4567
[19,] -7922.9295 -5046.1118
[20,] -6992.0713 -7922.9295
[21,] -4931.3682 -6992.0713
[22,] -7149.1548 -4931.3682
[23,] -4643.5113 -7149.1548
[24,] -1939.7481 -4643.5113
[25,] 225.3684 -1939.7481
[26,] 1510.0047 225.3684
[27,] 1735.9171 1510.0047
[28,] 2634.2944 1735.9171
[29,] 3166.0596 2634.2944
[30,] -3229.9164 3166.0596
[31,] -9358.7269 -3229.9164
[32,] -9271.3272 -9358.7269
[33,] -4957.2711 -9271.3272
[34,] -2659.0538 -4957.2711
[35,] 578.4296 -2659.0538
[36,] 3865.0859 578.4296
[37,] 6364.0747 3865.0859
[38,] 6418.3155 6364.0747
[39,] 9290.6561 6418.3155
[40,] 9304.5181 9290.6561
[41,] 9391.0709 9304.5181
[42,] 2208.3105 9391.0709
[43,] -1753.0569 2208.3105
[44,] -253.5178 -1753.0569
[45,] 2016.1708 -253.5178
[46,] 1445.9445 2016.1708
[47,] 4333.9427 1445.9445
[48,] 7607.3248 4333.9427
[49,] 8891.9427 7607.3248
[50,] 11817.2091 8891.9427
[51,] 9471.2429 11817.2091
[52,] 10732.2723 9471.2429
[53,] 10534.8281 10732.2723
[54,] 182.5911 10534.8281
[55,] -8097.3375 182.5911
[56,] -6286.4022 -8097.3375
[57,] 1818.6872 -6286.4022
[58,] 1428.2413 1818.6872
[59,] 5405.0039 1428.2413
[60,] 5061.8302 5405.0039
[61,] 6603.8378 5061.8302
[62,] 6922.8692 6603.8378
[63,] 4713.3807 6922.8692
[64,] 5814.6168 4713.3807
[65,] 4688.9749 5814.6168
[66,] -4421.2190 4688.9749
[67,] -9899.5372 -4421.2190
[68,] -10174.7977 -9899.5372
[69,] -10002.1157 -10174.7977
[70,] -11055.6971 -10002.1157
[71,] -7457.6259 -11055.6971
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 3945.8338 4965.9573
2 929.5462 3945.8338
3 1891.7464 929.5462
4 571.2728 1891.7464
5 860.3727 571.2728
6 -480.3222 860.3727
7 -2157.2873 -480.3222
8 -6529.9798 -2157.2873
9 -5680.3525 -6529.9798
10 -6563.5966 -5680.3525
11 -5030.5975 -6563.5966
12 -1743.8083 -5030.5975
13 -1760.9994 -1743.8083
14 -3926.2983 -1760.9994
15 -2769.8502 -3926.2983
16 -2164.6998 -2769.8502
17 -3037.4567 -2164.6998
18 -5046.1118 -3037.4567
19 -7922.9295 -5046.1118
20 -6992.0713 -7922.9295
21 -4931.3682 -6992.0713
22 -7149.1548 -4931.3682
23 -4643.5113 -7149.1548
24 -1939.7481 -4643.5113
25 225.3684 -1939.7481
26 1510.0047 225.3684
27 1735.9171 1510.0047
28 2634.2944 1735.9171
29 3166.0596 2634.2944
30 -3229.9164 3166.0596
31 -9358.7269 -3229.9164
32 -9271.3272 -9358.7269
33 -4957.2711 -9271.3272
34 -2659.0538 -4957.2711
35 578.4296 -2659.0538
36 3865.0859 578.4296
37 6364.0747 3865.0859
38 6418.3155 6364.0747
39 9290.6561 6418.3155
40 9304.5181 9290.6561
41 9391.0709 9304.5181
42 2208.3105 9391.0709
43 -1753.0569 2208.3105
44 -253.5178 -1753.0569
45 2016.1708 -253.5178
46 1445.9445 2016.1708
47 4333.9427 1445.9445
48 7607.3248 4333.9427
49 8891.9427 7607.3248
50 11817.2091 8891.9427
51 9471.2429 11817.2091
52 10732.2723 9471.2429
53 10534.8281 10732.2723
54 182.5911 10534.8281
55 -8097.3375 182.5911
56 -6286.4022 -8097.3375
57 1818.6872 -6286.4022
58 1428.2413 1818.6872
59 5405.0039 1428.2413
60 5061.8302 5405.0039
61 6603.8378 5061.8302
62 6922.8692 6603.8378
63 4713.3807 6922.8692
64 5814.6168 4713.3807
65 4688.9749 5814.6168
66 -4421.2190 4688.9749
67 -9899.5372 -4421.2190
68 -10174.7977 -9899.5372
69 -10002.1157 -10174.7977
70 -11055.6971 -10002.1157
71 -7457.6259 -11055.6971
> 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/7k5es1258484720.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/8glcc1258484720.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/9pohb1258484720.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/10bhz71258484720.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/11d0x81258484720.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/127qpa1258484720.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/134elp1258484720.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/14qosg1258484720.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/15mu3r1258484720.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/16gebn1258484720.tab")
+ }
>
> system("convert tmp/1zr841258484720.ps tmp/1zr841258484720.png")
> system("convert tmp/2lidu1258484720.ps tmp/2lidu1258484720.png")
> system("convert tmp/391qq1258484720.ps tmp/391qq1258484720.png")
> system("convert tmp/48xiu1258484720.ps tmp/48xiu1258484720.png")
> system("convert tmp/5pxk41258484720.ps tmp/5pxk41258484720.png")
> system("convert tmp/6fbjy1258484720.ps tmp/6fbjy1258484720.png")
> system("convert tmp/7k5es1258484720.ps tmp/7k5es1258484720.png")
> system("convert tmp/8glcc1258484720.ps tmp/8glcc1258484720.png")
> system("convert tmp/9pohb1258484720.ps tmp/9pohb1258484720.png")
> system("convert tmp/10bhz71258484720.ps tmp/10bhz71258484720.png")
>
>
> proc.time()
user system elapsed
2.641 1.649 6.301