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 = '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 t
1 562325 0 1 0 0 0 0 0 0 0 0 0 0 1
2 560854 0 0 1 0 0 0 0 0 0 0 0 0 2
3 555332 0 0 0 1 0 0 0 0 0 0 0 0 3
4 543599 0 0 0 0 1 0 0 0 0 0 0 0 4
5 536662 0 0 0 0 0 1 0 0 0 0 0 0 5
6 542722 0 0 0 0 0 0 1 0 0 0 0 0 6
7 593530 0 0 0 0 0 0 0 1 0 0 0 0 7
8 610763 0 0 0 0 0 0 0 0 1 0 0 0 8
9 612613 0 0 0 0 0 0 0 0 0 1 0 0 9
10 611324 0 0 0 0 0 0 0 0 0 0 1 0 10
11 594167 0 0 0 0 0 0 0 0 0 0 0 1 11
12 595454 0 0 0 0 0 0 0 0 0 0 0 0 12
13 590865 0 1 0 0 0 0 0 0 0 0 0 0 13
14 589379 0 0 1 0 0 0 0 0 0 0 0 0 14
15 584428 0 0 0 1 0 0 0 0 0 0 0 0 15
16 573100 0 0 0 0 1 0 0 0 0 0 0 0 16
17 567456 0 0 0 0 0 1 0 0 0 0 0 0 17
18 569028 0 0 0 0 0 0 1 0 0 0 0 0 18
19 620735 0 0 0 0 0 0 0 1 0 0 0 0 19
20 628884 0 0 0 0 0 0 0 0 1 0 0 0 20
21 628232 0 0 0 0 0 0 0 0 0 1 0 0 21
22 612117 0 0 0 0 0 0 0 0 0 0 1 0 22
23 595404 0 0 0 0 0 0 0 0 0 0 0 1 23
24 597141 0 0 0 0 0 0 0 0 0 0 0 0 24
25 593408 0 1 0 0 0 0 0 0 0 0 0 0 25
26 590072 0 0 1 0 0 0 0 0 0 0 0 0 26
27 579799 0 0 0 1 0 0 0 0 0 0 0 0 27
28 574205 0 0 0 0 1 0 0 0 0 0 0 0 28
29 572775 0 0 0 0 0 1 0 0 0 0 0 0 29
30 572942 0 0 0 0 0 0 1 0 0 0 0 0 30
31 619567 0 0 0 0 0 0 0 1 0 0 0 0 31
32 625809 0 0 0 0 0 0 0 0 1 0 0 0 32
33 619916 0 0 0 0 0 0 0 0 0 1 0 0 33
34 587625 0 0 0 0 0 0 0 0 0 0 1 0 34
35 565742 0 0 0 0 0 0 0 0 0 0 0 1 35
36 557274 0 0 0 0 0 0 0 0 0 0 0 0 36
37 560576 0 1 0 0 0 0 0 0 0 0 0 0 37
38 548854 0 0 1 0 0 0 0 0 0 0 0 0 38
39 531673 0 0 0 1 0 0 0 0 0 0 0 0 39
40 525919 0 0 0 0 1 0 0 0 0 0 0 0 40
41 511038 0 0 0 0 0 1 0 0 0 0 0 0 41
42 498662 1 0 0 0 0 0 1 0 0 0 0 0 42
43 555362 1 0 0 0 0 0 0 1 0 0 0 0 43
44 564591 1 0 0 0 0 0 0 0 1 0 0 0 44
45 541657 1 0 0 0 0 0 0 0 0 1 0 0 45
46 527070 1 0 0 0 0 0 0 0 0 0 1 0 46
47 509846 1 0 0 0 0 0 0 0 0 0 0 1 47
48 514258 1 0 0 0 0 0 0 0 0 0 0 0 48
49 516922 1 1 0 0 0 0 0 0 0 0 0 0 49
50 507561 1 0 1 0 0 0 0 0 0 0 0 0 50
51 492622 1 0 0 1 0 0 0 0 0 0 0 0 51
52 490243 1 0 0 0 1 0 0 0 0 0 0 0 52
53 469357 1 0 0 0 0 1 0 0 0 0 0 0 53
54 477580 1 0 0 0 0 0 1 0 0 0 0 0 54
55 528379 1 0 0 0 0 0 0 1 0 0 0 0 55
56 533590 1 0 0 0 0 0 0 0 1 0 0 0 56
57 517945 1 0 0 0 0 0 0 0 0 1 0 0 57
58 506174 1 0 0 0 0 0 0 0 0 0 1 0 58
59 501866 1 0 0 0 0 0 0 0 0 0 0 1 59
60 516141 1 0 0 0 0 0 0 0 0 0 0 0 60
61 528222 1 1 0 0 0 0 0 0 0 0 0 0 61
62 532638 1 0 1 0 0 0 0 0 0 0 0 0 62
63 536322 1 0 0 1 0 0 0 0 0 0 0 0 63
64 536535 1 0 0 0 1 0 0 0 0 0 0 0 64
65 523597 1 0 0 0 0 1 0 0 0 0 0 0 65
66 536214 1 0 0 0 0 0 1 0 0 0 0 0 66
67 586570 1 0 0 0 0 0 0 1 0 0 0 0 67
68 596594 1 0 0 0 0 0 0 0 1 0 0 0 68
69 580523 1 0 0 0 0 0 0 0 0 1 0 0 69
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
575726.7 -66470.4 -804.9 -4823.6 -13212.7 -19500.6
M5 M6 M7 M8 M9 M10
-30145.4 -16548.6 34425.2 43581.1 33498.2 13192.6
M11 t
-2456.5 192.1
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-42419 -13218 1659 15787 34486
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 575726.7 11286.3 51.011 < 2e-16 ***
X -66470.4 10150.7 -6.548 2.07e-08 ***
M1 -804.9 13039.6 -0.062 0.95101
M2 -4823.6 13029.1 -0.370 0.71264
M3 -13212.7 13023.5 -1.015 0.31477
M4 -19500.6 13022.6 -1.497 0.13999
M5 -30145.4 13026.5 -2.314 0.02442 *
M6 -16548.6 13057.0 -1.267 0.21035
M7 34425.2 13042.9 2.639 0.01078 *
M8 43581.1 13033.5 3.344 0.00149 **
M9 33498.2 13028.9 2.571 0.01287 *
M10 13192.6 13605.5 0.970 0.33647
M11 -2456.5 13598.6 -0.181 0.85731
t 192.1 249.2 0.771 0.44403
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21500 on 55 degrees of freedom
Multiple R-squared: 0.7713, Adjusted R-squared: 0.7172
F-statistic: 14.27 on 13 and 55 DF, p-value: 3.21e-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,] 4.753833e-05 9.507666e-05 0.99995246
[2,] 1.197065e-05 2.394130e-05 0.99998803
[3,] 7.672161e-07 1.534432e-06 0.99999923
[4,] 1.720375e-05 3.440751e-05 0.99998280
[5,] 2.547572e-05 5.095144e-05 0.99997452
[6,] 5.228037e-04 1.045607e-03 0.99947720
[7,] 1.007518e-03 2.015036e-03 0.99899248
[8,] 1.169305e-03 2.338610e-03 0.99883070
[9,] 7.505931e-04 1.501186e-03 0.99924941
[10,] 5.030646e-04 1.006129e-03 0.99949694
[11,] 4.568671e-04 9.137342e-04 0.99954313
[12,] 2.574670e-04 5.149340e-04 0.99974253
[13,] 2.184938e-04 4.369876e-04 0.99978151
[14,] 1.155182e-04 2.310364e-04 0.99988448
[15,] 7.085767e-05 1.417153e-04 0.99992914
[16,] 7.177909e-05 1.435582e-04 0.99992822
[17,] 2.319223e-04 4.638446e-04 0.99976808
[18,] 6.139203e-03 1.227841e-02 0.99386080
[19,] 3.750524e-02 7.501047e-02 0.96249476
[20,] 1.197016e-01 2.394032e-01 0.88029839
[21,] 1.440860e-01 2.881721e-01 0.85591395
[22,] 1.798706e-01 3.597413e-01 0.82012937
[23,] 2.279711e-01 4.559421e-01 0.77202893
[24,] 2.295511e-01 4.591021e-01 0.77044893
[25,] 2.481994e-01 4.963988e-01 0.75180061
[26,] 2.006970e-01 4.013941e-01 0.79930296
[27,] 1.849683e-01 3.699366e-01 0.81503169
[28,] 2.012270e-01 4.024540e-01 0.79877298
[29,] 2.436471e-01 4.872941e-01 0.75635294
[30,] 4.229671e-01 8.459342e-01 0.57703290
[31,] 6.096268e-01 7.807464e-01 0.39037322
[32,] 7.916930e-01 4.166141e-01 0.20830703
[33,] 9.374240e-01 1.251521e-01 0.06257605
[34,] 9.895938e-01 2.081234e-02 0.01040617
[35,] 9.886172e-01 2.276569e-02 0.01138284
[36,] 9.934168e-01 1.316646e-02 0.00658323
> postscript(file="/var/www/html/rcomp/tmp/1wthc1258711785.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/2rp9d1258711785.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/378vs1258711785.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/41k491258711785.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/5oi281258711785.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
-12788.9560 -10433.2893 -7758.2893 -13395.4560 -9879.7893 -17608.6931
7 8 9 10 11 12
-17966.5265 -10081.5265 1659.3069 20483.8455 18783.8455 17422.2455
13 14 15 16 17 18
13446.0367 15786.7034 19032.7034 13800.5367 18609.2034 6392.2996
19 20 21 22 23 24
6933.4663 5734.4663 14973.2996 18971.8382 17715.8382 16804.2382
25 26 27 28 29 30
13684.0294 14174.6961 12098.6961 12600.5294 21623.1961 8001.2923
31 32 33 34 35 36
3460.4590 354.4590 4352.2923 -7825.1691 -14251.1691 -25367.7691
37 38 39 40 41 42
-21452.9779 -29348.3112 -38332.3112 -37990.4779 -42418.8112 -2113.2923
43 44 45 46 47 48
3420.8744 3301.8744 -9741.2923 -4214.7537 -5981.7537 -4218.3537
49 50 51 52 53 54
-941.5625 -6475.8958 -13217.8958 -9501.0625 -19934.3958 -25500.2996
55 56 57 58 59 60
-25867.1329 -30004.1329 -35758.2996 -27415.7610 -16266.7610 -4640.3610
61 62 63 64 65 66
8053.4302 16296.0969 28177.0969 34485.9302 32000.5969 30828.6931
67 68 69
30018.8598 30694.8598 24514.6931
> postscript(file="/var/www/html/rcomp/tmp/69s3b1258711785.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 -12788.9560 NA
1 -10433.2893 -12788.9560
2 -7758.2893 -10433.2893
3 -13395.4560 -7758.2893
4 -9879.7893 -13395.4560
5 -17608.6931 -9879.7893
6 -17966.5265 -17608.6931
7 -10081.5265 -17966.5265
8 1659.3069 -10081.5265
9 20483.8455 1659.3069
10 18783.8455 20483.8455
11 17422.2455 18783.8455
12 13446.0367 17422.2455
13 15786.7034 13446.0367
14 19032.7034 15786.7034
15 13800.5367 19032.7034
16 18609.2034 13800.5367
17 6392.2996 18609.2034
18 6933.4663 6392.2996
19 5734.4663 6933.4663
20 14973.2996 5734.4663
21 18971.8382 14973.2996
22 17715.8382 18971.8382
23 16804.2382 17715.8382
24 13684.0294 16804.2382
25 14174.6961 13684.0294
26 12098.6961 14174.6961
27 12600.5294 12098.6961
28 21623.1961 12600.5294
29 8001.2923 21623.1961
30 3460.4590 8001.2923
31 354.4590 3460.4590
32 4352.2923 354.4590
33 -7825.1691 4352.2923
34 -14251.1691 -7825.1691
35 -25367.7691 -14251.1691
36 -21452.9779 -25367.7691
37 -29348.3112 -21452.9779
38 -38332.3112 -29348.3112
39 -37990.4779 -38332.3112
40 -42418.8112 -37990.4779
41 -2113.2923 -42418.8112
42 3420.8744 -2113.2923
43 3301.8744 3420.8744
44 -9741.2923 3301.8744
45 -4214.7537 -9741.2923
46 -5981.7537 -4214.7537
47 -4218.3537 -5981.7537
48 -941.5625 -4218.3537
49 -6475.8958 -941.5625
50 -13217.8958 -6475.8958
51 -9501.0625 -13217.8958
52 -19934.3958 -9501.0625
53 -25500.2996 -19934.3958
54 -25867.1329 -25500.2996
55 -30004.1329 -25867.1329
56 -35758.2996 -30004.1329
57 -27415.7610 -35758.2996
58 -16266.7610 -27415.7610
59 -4640.3610 -16266.7610
60 8053.4302 -4640.3610
61 16296.0969 8053.4302
62 28177.0969 16296.0969
63 34485.9302 28177.0969
64 32000.5969 34485.9302
65 30828.6931 32000.5969
66 30018.8598 30828.6931
67 30694.8598 30018.8598
68 24514.6931 30694.8598
69 NA 24514.6931
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -10433.2893 -12788.9560
[2,] -7758.2893 -10433.2893
[3,] -13395.4560 -7758.2893
[4,] -9879.7893 -13395.4560
[5,] -17608.6931 -9879.7893
[6,] -17966.5265 -17608.6931
[7,] -10081.5265 -17966.5265
[8,] 1659.3069 -10081.5265
[9,] 20483.8455 1659.3069
[10,] 18783.8455 20483.8455
[11,] 17422.2455 18783.8455
[12,] 13446.0367 17422.2455
[13,] 15786.7034 13446.0367
[14,] 19032.7034 15786.7034
[15,] 13800.5367 19032.7034
[16,] 18609.2034 13800.5367
[17,] 6392.2996 18609.2034
[18,] 6933.4663 6392.2996
[19,] 5734.4663 6933.4663
[20,] 14973.2996 5734.4663
[21,] 18971.8382 14973.2996
[22,] 17715.8382 18971.8382
[23,] 16804.2382 17715.8382
[24,] 13684.0294 16804.2382
[25,] 14174.6961 13684.0294
[26,] 12098.6961 14174.6961
[27,] 12600.5294 12098.6961
[28,] 21623.1961 12600.5294
[29,] 8001.2923 21623.1961
[30,] 3460.4590 8001.2923
[31,] 354.4590 3460.4590
[32,] 4352.2923 354.4590
[33,] -7825.1691 4352.2923
[34,] -14251.1691 -7825.1691
[35,] -25367.7691 -14251.1691
[36,] -21452.9779 -25367.7691
[37,] -29348.3112 -21452.9779
[38,] -38332.3112 -29348.3112
[39,] -37990.4779 -38332.3112
[40,] -42418.8112 -37990.4779
[41,] -2113.2923 -42418.8112
[42,] 3420.8744 -2113.2923
[43,] 3301.8744 3420.8744
[44,] -9741.2923 3301.8744
[45,] -4214.7537 -9741.2923
[46,] -5981.7537 -4214.7537
[47,] -4218.3537 -5981.7537
[48,] -941.5625 -4218.3537
[49,] -6475.8958 -941.5625
[50,] -13217.8958 -6475.8958
[51,] -9501.0625 -13217.8958
[52,] -19934.3958 -9501.0625
[53,] -25500.2996 -19934.3958
[54,] -25867.1329 -25500.2996
[55,] -30004.1329 -25867.1329
[56,] -35758.2996 -30004.1329
[57,] -27415.7610 -35758.2996
[58,] -16266.7610 -27415.7610
[59,] -4640.3610 -16266.7610
[60,] 8053.4302 -4640.3610
[61,] 16296.0969 8053.4302
[62,] 28177.0969 16296.0969
[63,] 34485.9302 28177.0969
[64,] 32000.5969 34485.9302
[65,] 30828.6931 32000.5969
[66,] 30018.8598 30828.6931
[67,] 30694.8598 30018.8598
[68,] 24514.6931 30694.8598
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -10433.2893 -12788.9560
2 -7758.2893 -10433.2893
3 -13395.4560 -7758.2893
4 -9879.7893 -13395.4560
5 -17608.6931 -9879.7893
6 -17966.5265 -17608.6931
7 -10081.5265 -17966.5265
8 1659.3069 -10081.5265
9 20483.8455 1659.3069
10 18783.8455 20483.8455
11 17422.2455 18783.8455
12 13446.0367 17422.2455
13 15786.7034 13446.0367
14 19032.7034 15786.7034
15 13800.5367 19032.7034
16 18609.2034 13800.5367
17 6392.2996 18609.2034
18 6933.4663 6392.2996
19 5734.4663 6933.4663
20 14973.2996 5734.4663
21 18971.8382 14973.2996
22 17715.8382 18971.8382
23 16804.2382 17715.8382
24 13684.0294 16804.2382
25 14174.6961 13684.0294
26 12098.6961 14174.6961
27 12600.5294 12098.6961
28 21623.1961 12600.5294
29 8001.2923 21623.1961
30 3460.4590 8001.2923
31 354.4590 3460.4590
32 4352.2923 354.4590
33 -7825.1691 4352.2923
34 -14251.1691 -7825.1691
35 -25367.7691 -14251.1691
36 -21452.9779 -25367.7691
37 -29348.3112 -21452.9779
38 -38332.3112 -29348.3112
39 -37990.4779 -38332.3112
40 -42418.8112 -37990.4779
41 -2113.2923 -42418.8112
42 3420.8744 -2113.2923
43 3301.8744 3420.8744
44 -9741.2923 3301.8744
45 -4214.7537 -9741.2923
46 -5981.7537 -4214.7537
47 -4218.3537 -5981.7537
48 -941.5625 -4218.3537
49 -6475.8958 -941.5625
50 -13217.8958 -6475.8958
51 -9501.0625 -13217.8958
52 -19934.3958 -9501.0625
53 -25500.2996 -19934.3958
54 -25867.1329 -25500.2996
55 -30004.1329 -25867.1329
56 -35758.2996 -30004.1329
57 -27415.7610 -35758.2996
58 -16266.7610 -27415.7610
59 -4640.3610 -16266.7610
60 8053.4302 -4640.3610
61 16296.0969 8053.4302
62 28177.0969 16296.0969
63 34485.9302 28177.0969
64 32000.5969 34485.9302
65 30828.6931 32000.5969
66 30018.8598 30828.6931
67 30694.8598 30018.8598
68 24514.6931 30694.8598
> 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/7c0jc1258711785.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/8qnez1258711785.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/9e7ha1258711785.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/10jgr41258711785.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/11rd571258711785.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/12eiqy1258711785.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/13e41w1258711786.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/14lgcr1258711786.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/15p4dp1258711786.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/16asy81258711786.tab")
+ }
> system("convert tmp/1wthc1258711785.ps tmp/1wthc1258711785.png")
> system("convert tmp/2rp9d1258711785.ps tmp/2rp9d1258711785.png")
> system("convert tmp/378vs1258711785.ps tmp/378vs1258711785.png")
> system("convert tmp/41k491258711785.ps tmp/41k491258711785.png")
> system("convert tmp/5oi281258711785.ps tmp/5oi281258711785.png")
> system("convert tmp/69s3b1258711785.ps tmp/69s3b1258711785.png")
> system("convert tmp/7c0jc1258711785.ps tmp/7c0jc1258711785.png")
> system("convert tmp/8qnez1258711785.ps tmp/8qnez1258711785.png")
> system("convert tmp/9e7ha1258711785.ps tmp/9e7ha1258711785.png")
> system("convert tmp/10jgr41258711785.ps tmp/10jgr41258711785.png")
>
>
> proc.time()
user system elapsed
2.477 1.550 2.952