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(562325,0,560854,0,555332,0,543599,0,536662,0,542722,0,593530,0,610763,0,612613,0,611324,0,594167,0,595454,0,590865,0,589379,0,584428,0,573100,0,567456,0,569028,0,620735,0,628884,0,628232,0,612117,0,595404,0,597141,0,593408,0,590072,0,579799,0,574205,0,572775,0,572942,0,619567,0,625809,0,619916,0,587625,0,565742,0,557274,0,560576,0,548854,0,531673,0,525919,0,511038,0,498662,1,555362,1,564591,1,541657,1,527070,1,509846,1,514258,1,516922,1,507561,1,492622,1,490243,1,469357,1,477580,1,528379,1,533590,1,517945,1,506174,1,501866,1,516141,1,528222,1,532638,1,536322,1,536535,1,523597,1,536214,1,586570,1,596594,1,580523,1),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 562325 0 1 0 0 0 0 0 0 0 0 0 0
2 560854 0 0 1 0 0 0 0 0 0 0 0 0
3 555332 0 0 0 1 0 0 0 0 0 0 0 0
4 543599 0 0 0 0 1 0 0 0 0 0 0 0
5 536662 0 0 0 0 0 1 0 0 0 0 0 0
6 542722 0 0 0 0 0 0 1 0 0 0 0 0
7 593530 0 0 0 0 0 0 0 1 0 0 0 0
8 610763 0 0 0 0 0 0 0 0 1 0 0 0
9 612613 0 0 0 0 0 0 0 0 0 1 0 0
10 611324 0 0 0 0 0 0 0 0 0 0 1 0
11 594167 0 0 0 0 0 0 0 0 0 0 0 1
12 595454 0 0 0 0 0 0 0 0 0 0 0 0
13 590865 0 1 0 0 0 0 0 0 0 0 0 0
14 589379 0 0 1 0 0 0 0 0 0 0 0 0
15 584428 0 0 0 1 0 0 0 0 0 0 0 0
16 573100 0 0 0 0 1 0 0 0 0 0 0 0
17 567456 0 0 0 0 0 1 0 0 0 0 0 0
18 569028 0 0 0 0 0 0 1 0 0 0 0 0
19 620735 0 0 0 0 0 0 0 1 0 0 0 0
20 628884 0 0 0 0 0 0 0 0 1 0 0 0
21 628232 0 0 0 0 0 0 0 0 0 1 0 0
22 612117 0 0 0 0 0 0 0 0 0 0 1 0
23 595404 0 0 0 0 0 0 0 0 0 0 0 1
24 597141 0 0 0 0 0 0 0 0 0 0 0 0
25 593408 0 1 0 0 0 0 0 0 0 0 0 0
26 590072 0 0 1 0 0 0 0 0 0 0 0 0
27 579799 0 0 0 1 0 0 0 0 0 0 0 0
28 574205 0 0 0 0 1 0 0 0 0 0 0 0
29 572775 0 0 0 0 0 1 0 0 0 0 0 0
30 572942 0 0 0 0 0 0 1 0 0 0 0 0
31 619567 0 0 0 0 0 0 0 1 0 0 0 0
32 625809 0 0 0 0 0 0 0 0 1 0 0 0
33 619916 0 0 0 0 0 0 0 0 0 1 0 0
34 587625 0 0 0 0 0 0 0 0 0 0 1 0
35 565742 0 0 0 0 0 0 0 0 0 0 0 1
36 557274 0 0 0 0 0 0 0 0 0 0 0 0
37 560576 0 1 0 0 0 0 0 0 0 0 0 0
38 548854 0 0 1 0 0 0 0 0 0 0 0 0
39 531673 0 0 0 1 0 0 0 0 0 0 0 0
40 525919 0 0 0 0 1 0 0 0 0 0 0 0
41 511038 0 0 0 0 0 1 0 0 0 0 0 0
42 498662 1 0 0 0 0 0 1 0 0 0 0 0
43 555362 1 0 0 0 0 0 0 1 0 0 0 0
44 564591 1 0 0 0 0 0 0 0 1 0 0 0
45 541657 1 0 0 0 0 0 0 0 0 1 0 0
46 527070 1 0 0 0 0 0 0 0 0 0 1 0
47 509846 1 0 0 0 0 0 0 0 0 0 0 1
48 514258 1 0 0 0 0 0 0 0 0 0 0 0
49 516922 1 1 0 0 0 0 0 0 0 0 0 0
50 507561 1 0 1 0 0 0 0 0 0 0 0 0
51 492622 1 0 0 1 0 0 0 0 0 0 0 0
52 490243 1 0 0 0 1 0 0 0 0 0 0 0
53 469357 1 0 0 0 0 1 0 0 0 0 0 0
54 477580 1 0 0 0 0 0 1 0 0 0 0 0
55 528379 1 0 0 0 0 0 0 1 0 0 0 0
56 533590 1 0 0 0 0 0 0 0 1 0 0 0
57 517945 1 0 0 0 0 0 0 0 0 1 0 0
58 506174 1 0 0 0 0 0 0 0 0 0 1 0
59 501866 1 0 0 0 0 0 0 0 0 0 0 1
60 516141 1 0 0 0 0 0 0 0 0 0 0 0
61 528222 1 1 0 0 0 0 0 0 0 0 0 0
62 532638 1 0 1 0 0 0 0 0 0 0 0 0
63 536322 1 0 0 1 0 0 0 0 0 0 0 0
64 536535 1 0 0 0 1 0 0 0 0 0 0 0
65 523597 1 0 0 0 0 1 0 0 0 0 0 0
66 536214 1 0 0 0 0 0 1 0 0 0 0 0
67 586570 1 0 0 0 0 0 0 1 0 0 0 0
68 596594 1 0 0 0 0 0 0 0 1 0 0 0
69 580523 1 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
579978 -59810 -1321 -5148 -13345 -19441
M5 M6 M7 M8 M9 M10
-29893 -17215 33951 43299 33408 12808
M11
-2649
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-39046 -14200 1124 15242 35808
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 579978 9812 59.110 < 2e-16 ***
X -59810 5311 -11.262 5.18e-16 ***
M1 -1321 12975 -0.102 0.91925
M2 -5148 12975 -0.397 0.69305
M3 -13345 12975 -1.029 0.30813
M4 -19441 12975 -1.498 0.13967
M5 -29894 12975 -2.304 0.02496 *
M6 -17215 12981 -1.326 0.19018
M7 33951 12981 2.615 0.01143 *
M8 43299 12981 3.336 0.00152 **
M9 33408 12981 2.574 0.01274 *
M10 12808 13547 0.945 0.34848
M11 -2649 13547 -0.196 0.84570
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21420 on 56 degrees of freedom
Multiple R-squared: 0.7688, Adjusted R-squared: 0.7193
F-statistic: 15.52 on 12 and 56 DF, p-value: 1.043e-13
> 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.6139351626 0.7721296748 0.3860648
[2,] 0.5736945609 0.8526108783 0.4263054
[3,] 0.5044652212 0.9910695576 0.4955348
[4,] 0.4473676287 0.8947352575 0.5526324
[5,] 0.3558179475 0.7116358950 0.6441821
[6,] 0.2766412760 0.5532825520 0.7233587
[7,] 0.2058454088 0.4116908177 0.7941546
[8,] 0.1495006886 0.2990013773 0.8504993
[9,] 0.1065109988 0.2130219977 0.8934890
[10,] 0.0850573626 0.1701147251 0.9149426
[11,] 0.0666413727 0.1332827453 0.9333586
[12,] 0.0478024954 0.0956049908 0.9521975
[13,] 0.0376025644 0.0752051289 0.9623974
[14,] 0.0415001544 0.0830003087 0.9584998
[15,] 0.0341952294 0.0683904589 0.9658048
[16,] 0.0245122751 0.0490245503 0.9754877
[17,] 0.0158655467 0.0317310934 0.9841345
[18,] 0.0125958449 0.0251916899 0.9874042
[19,] 0.0167071761 0.0334143522 0.9832928
[20,] 0.0234887796 0.0469775592 0.9765112
[21,] 0.0378405963 0.0756811925 0.9621594
[22,] 0.0332138015 0.0664276029 0.9667862
[23,] 0.0358881320 0.0717762639 0.9641119
[24,] 0.0491743545 0.0983487089 0.9508256
[25,] 0.0536812492 0.1073624985 0.9463188
[26,] 0.0716896034 0.1433792068 0.9283104
[27,] 0.0463694012 0.0927388024 0.9536306
[28,] 0.0283293212 0.0566586424 0.9716707
[29,] 0.0162909670 0.0325819340 0.9837090
[30,] 0.0101245442 0.0202490884 0.9898755
[31,] 0.0062248839 0.0124497679 0.9937751
[32,] 0.0031892803 0.0063785606 0.9968107
[33,] 0.0014547471 0.0029094943 0.9985453
[34,] 0.0006557838 0.0013115676 0.9993442
[35,] 0.0003163843 0.0006327686 0.9996836
[36,] 0.0002342002 0.0004684004 0.9997658
[37,] 0.0001885569 0.0003771138 0.9998114
[38,] 0.0002325564 0.0004651128 0.9997674
> postscript(file="/var/www/html/rcomp/tmp/1deaq1258663180.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/200241258663180.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/3kr6p1258663180.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/45w3m1258663180.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/5dzka1258663180.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
-16331.4877 -13975.8210 -11300.8210 -16937.9877 -13422.3210 -20041.2316
7 8 9 10 11 12
-20399.0649 -12514.0649 -773.2316 18537.8148 16837.8148 15476.2148
13 14 15 16 17 18
12208.5123 14549.1790 17795.1790 12563.0123 17371.6790 6264.7684
19 20 21 22 23 24
6805.9351 5606.9351 14845.7684 19330.8148 18074.8148 17163.2148
25 26 27 28 29 30
14751.5123 15242.1790 13166.1790 13668.0123 22690.6790 10178.7684
31 32 33 34 35 36
5637.9351 2531.9351 6529.7684 -5161.1852 -11587.1852 -22703.7852
37 38 39 40 41 42
-18080.4877 -25975.8210 -34959.8210 -34617.9877 -39046.3210 -4290.7684
43 44 45 46 47 48
1243.3982 1124.3982 -11918.7684 -5905.7221 -7672.7221 -5909.3221
49 50 51 52 53 54
-1924.0246 -7458.3579 -14200.3579 -10483.5246 -20916.8579 -25372.7684
55 56 57 58 59 60
-25739.6018 -29876.6018 -35630.7684 -26801.7221 -15652.7221 -4026.3221
61 62 63 64 65 66
9375.9754 17618.6421 29499.6421 35808.4754 33323.1421 33261.2316
67 68 69
32451.3982 33127.3982 26947.2316
> postscript(file="/var/www/html/rcomp/tmp/603el1258663180.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 -16331.4877 NA
1 -13975.8210 -16331.4877
2 -11300.8210 -13975.8210
3 -16937.9877 -11300.8210
4 -13422.3210 -16937.9877
5 -20041.2316 -13422.3210
6 -20399.0649 -20041.2316
7 -12514.0649 -20399.0649
8 -773.2316 -12514.0649
9 18537.8148 -773.2316
10 16837.8148 18537.8148
11 15476.2148 16837.8148
12 12208.5123 15476.2148
13 14549.1790 12208.5123
14 17795.1790 14549.1790
15 12563.0123 17795.1790
16 17371.6790 12563.0123
17 6264.7684 17371.6790
18 6805.9351 6264.7684
19 5606.9351 6805.9351
20 14845.7684 5606.9351
21 19330.8148 14845.7684
22 18074.8148 19330.8148
23 17163.2148 18074.8148
24 14751.5123 17163.2148
25 15242.1790 14751.5123
26 13166.1790 15242.1790
27 13668.0123 13166.1790
28 22690.6790 13668.0123
29 10178.7684 22690.6790
30 5637.9351 10178.7684
31 2531.9351 5637.9351
32 6529.7684 2531.9351
33 -5161.1852 6529.7684
34 -11587.1852 -5161.1852
35 -22703.7852 -11587.1852
36 -18080.4877 -22703.7852
37 -25975.8210 -18080.4877
38 -34959.8210 -25975.8210
39 -34617.9877 -34959.8210
40 -39046.3210 -34617.9877
41 -4290.7684 -39046.3210
42 1243.3982 -4290.7684
43 1124.3982 1243.3982
44 -11918.7684 1124.3982
45 -5905.7221 -11918.7684
46 -7672.7221 -5905.7221
47 -5909.3221 -7672.7221
48 -1924.0246 -5909.3221
49 -7458.3579 -1924.0246
50 -14200.3579 -7458.3579
51 -10483.5246 -14200.3579
52 -20916.8579 -10483.5246
53 -25372.7684 -20916.8579
54 -25739.6018 -25372.7684
55 -29876.6018 -25739.6018
56 -35630.7684 -29876.6018
57 -26801.7221 -35630.7684
58 -15652.7221 -26801.7221
59 -4026.3221 -15652.7221
60 9375.9754 -4026.3221
61 17618.6421 9375.9754
62 29499.6421 17618.6421
63 35808.4754 29499.6421
64 33323.1421 35808.4754
65 33261.2316 33323.1421
66 32451.3982 33261.2316
67 33127.3982 32451.3982
68 26947.2316 33127.3982
69 NA 26947.2316
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -13975.8210 -16331.4877
[2,] -11300.8210 -13975.8210
[3,] -16937.9877 -11300.8210
[4,] -13422.3210 -16937.9877
[5,] -20041.2316 -13422.3210
[6,] -20399.0649 -20041.2316
[7,] -12514.0649 -20399.0649
[8,] -773.2316 -12514.0649
[9,] 18537.8148 -773.2316
[10,] 16837.8148 18537.8148
[11,] 15476.2148 16837.8148
[12,] 12208.5123 15476.2148
[13,] 14549.1790 12208.5123
[14,] 17795.1790 14549.1790
[15,] 12563.0123 17795.1790
[16,] 17371.6790 12563.0123
[17,] 6264.7684 17371.6790
[18,] 6805.9351 6264.7684
[19,] 5606.9351 6805.9351
[20,] 14845.7684 5606.9351
[21,] 19330.8148 14845.7684
[22,] 18074.8148 19330.8148
[23,] 17163.2148 18074.8148
[24,] 14751.5123 17163.2148
[25,] 15242.1790 14751.5123
[26,] 13166.1790 15242.1790
[27,] 13668.0123 13166.1790
[28,] 22690.6790 13668.0123
[29,] 10178.7684 22690.6790
[30,] 5637.9351 10178.7684
[31,] 2531.9351 5637.9351
[32,] 6529.7684 2531.9351
[33,] -5161.1852 6529.7684
[34,] -11587.1852 -5161.1852
[35,] -22703.7852 -11587.1852
[36,] -18080.4877 -22703.7852
[37,] -25975.8210 -18080.4877
[38,] -34959.8210 -25975.8210
[39,] -34617.9877 -34959.8210
[40,] -39046.3210 -34617.9877
[41,] -4290.7684 -39046.3210
[42,] 1243.3982 -4290.7684
[43,] 1124.3982 1243.3982
[44,] -11918.7684 1124.3982
[45,] -5905.7221 -11918.7684
[46,] -7672.7221 -5905.7221
[47,] -5909.3221 -7672.7221
[48,] -1924.0246 -5909.3221
[49,] -7458.3579 -1924.0246
[50,] -14200.3579 -7458.3579
[51,] -10483.5246 -14200.3579
[52,] -20916.8579 -10483.5246
[53,] -25372.7684 -20916.8579
[54,] -25739.6018 -25372.7684
[55,] -29876.6018 -25739.6018
[56,] -35630.7684 -29876.6018
[57,] -26801.7221 -35630.7684
[58,] -15652.7221 -26801.7221
[59,] -4026.3221 -15652.7221
[60,] 9375.9754 -4026.3221
[61,] 17618.6421 9375.9754
[62,] 29499.6421 17618.6421
[63,] 35808.4754 29499.6421
[64,] 33323.1421 35808.4754
[65,] 33261.2316 33323.1421
[66,] 32451.3982 33261.2316
[67,] 33127.3982 32451.3982
[68,] 26947.2316 33127.3982
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -13975.8210 -16331.4877
2 -11300.8210 -13975.8210
3 -16937.9877 -11300.8210
4 -13422.3210 -16937.9877
5 -20041.2316 -13422.3210
6 -20399.0649 -20041.2316
7 -12514.0649 -20399.0649
8 -773.2316 -12514.0649
9 18537.8148 -773.2316
10 16837.8148 18537.8148
11 15476.2148 16837.8148
12 12208.5123 15476.2148
13 14549.1790 12208.5123
14 17795.1790 14549.1790
15 12563.0123 17795.1790
16 17371.6790 12563.0123
17 6264.7684 17371.6790
18 6805.9351 6264.7684
19 5606.9351 6805.9351
20 14845.7684 5606.9351
21 19330.8148 14845.7684
22 18074.8148 19330.8148
23 17163.2148 18074.8148
24 14751.5123 17163.2148
25 15242.1790 14751.5123
26 13166.1790 15242.1790
27 13668.0123 13166.1790
28 22690.6790 13668.0123
29 10178.7684 22690.6790
30 5637.9351 10178.7684
31 2531.9351 5637.9351
32 6529.7684 2531.9351
33 -5161.1852 6529.7684
34 -11587.1852 -5161.1852
35 -22703.7852 -11587.1852
36 -18080.4877 -22703.7852
37 -25975.8210 -18080.4877
38 -34959.8210 -25975.8210
39 -34617.9877 -34959.8210
40 -39046.3210 -34617.9877
41 -4290.7684 -39046.3210
42 1243.3982 -4290.7684
43 1124.3982 1243.3982
44 -11918.7684 1124.3982
45 -5905.7221 -11918.7684
46 -7672.7221 -5905.7221
47 -5909.3221 -7672.7221
48 -1924.0246 -5909.3221
49 -7458.3579 -1924.0246
50 -14200.3579 -7458.3579
51 -10483.5246 -14200.3579
52 -20916.8579 -10483.5246
53 -25372.7684 -20916.8579
54 -25739.6018 -25372.7684
55 -29876.6018 -25739.6018
56 -35630.7684 -29876.6018
57 -26801.7221 -35630.7684
58 -15652.7221 -26801.7221
59 -4026.3221 -15652.7221
60 9375.9754 -4026.3221
61 17618.6421 9375.9754
62 29499.6421 17618.6421
63 35808.4754 29499.6421
64 33323.1421 35808.4754
65 33261.2316 33323.1421
66 32451.3982 33261.2316
67 33127.3982 32451.3982
68 26947.2316 33127.3982
> 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/7bnji1258663180.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/88sjp1258663180.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/95m821258663180.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/10dwll1258663180.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/11281v1258663180.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/1209zp1258663180.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/133lk81258663180.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/14kti11258663180.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/153q981258663180.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/16xajy1258663180.tab")
+ }
> system("convert tmp/1deaq1258663180.ps tmp/1deaq1258663180.png")
> system("convert tmp/200241258663180.ps tmp/200241258663180.png")
> system("convert tmp/3kr6p1258663180.ps tmp/3kr6p1258663180.png")
> system("convert tmp/45w3m1258663180.ps tmp/45w3m1258663180.png")
> system("convert tmp/5dzka1258663180.ps tmp/5dzka1258663180.png")
> system("convert tmp/603el1258663180.ps tmp/603el1258663180.png")
> system("convert tmp/7bnji1258663180.ps tmp/7bnji1258663180.png")
> system("convert tmp/88sjp1258663180.ps tmp/88sjp1258663180.png")
> system("convert tmp/95m821258663180.ps tmp/95m821258663180.png")
> system("convert tmp/10dwll1258663180.ps tmp/10dwll1258663180.png")
>
>
> proc.time()
user system elapsed
2.494 1.547 2.882