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(2849.27,10872,2921.44,10625,2981.85,10407,3080.58,10463,3106.22,10556,3119.31,10646,3061.26,10702,3097.31,11353,3161.69,11346,3257.16,11451,3277.01,11964,3295.32,12574,3363.99,13031,3494.17,13812,3667.03,14544,3813.06,14931,3917.96,14886,3895.51,16005,3801.06,17064,3570.12,15168,3701.61,16050,3862.27,15839,3970.1,15137,4138.52,14954,4199.75,15648,4290.89,15305,4443.91,15579,4502.64,16348,4356.98,15928,4591.27,16171,4696.96,15937,4621.4,15713,4562.84,15594,4202.52,15683,4296.49,16438,4435.23,17032,4105.18,17696,4116.68,17745,3844.49,19394,3720.98,20148,3674.4,20108,3857.62,18584,3801.06,18441,3504.37,18391,3032.6,19178,3047.03,18079,2962.34,18483,2197.82,19644,2014.45,19195,1862.83,19650,1905.41,20830,1810.99,23595,1670.07,22937,1864.44,21814,2052.02,21928,2029.6,21777,2070.83,21383,2293.41,21467,2443.27,22052,2513.17,22680),dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> 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 2849.27 10872 1 0 0 0 0 0 0 0 0 0 0
2 2921.44 10625 0 1 0 0 0 0 0 0 0 0 0
3 2981.85 10407 0 0 1 0 0 0 0 0 0 0 0
4 3080.58 10463 0 0 0 1 0 0 0 0 0 0 0
5 3106.22 10556 0 0 0 0 1 0 0 0 0 0 0
6 3119.31 10646 0 0 0 0 0 1 0 0 0 0 0
7 3061.26 10702 0 0 0 0 0 0 1 0 0 0 0
8 3097.31 11353 0 0 0 0 0 0 0 1 0 0 0
9 3161.69 11346 0 0 0 0 0 0 0 0 1 0 0
10 3257.16 11451 0 0 0 0 0 0 0 0 0 1 0
11 3277.01 11964 0 0 0 0 0 0 0 0 0 0 1
12 3295.32 12574 0 0 0 0 0 0 0 0 0 0 0
13 3363.99 13031 1 0 0 0 0 0 0 0 0 0 0
14 3494.17 13812 0 1 0 0 0 0 0 0 0 0 0
15 3667.03 14544 0 0 1 0 0 0 0 0 0 0 0
16 3813.06 14931 0 0 0 1 0 0 0 0 0 0 0
17 3917.96 14886 0 0 0 0 1 0 0 0 0 0 0
18 3895.51 16005 0 0 0 0 0 1 0 0 0 0 0
19 3801.06 17064 0 0 0 0 0 0 1 0 0 0 0
20 3570.12 15168 0 0 0 0 0 0 0 1 0 0 0
21 3701.61 16050 0 0 0 0 0 0 0 0 1 0 0
22 3862.27 15839 0 0 0 0 0 0 0 0 0 1 0
23 3970.10 15137 0 0 0 0 0 0 0 0 0 0 1
24 4138.52 14954 0 0 0 0 0 0 0 0 0 0 0
25 4199.75 15648 1 0 0 0 0 0 0 0 0 0 0
26 4290.89 15305 0 1 0 0 0 0 0 0 0 0 0
27 4443.91 15579 0 0 1 0 0 0 0 0 0 0 0
28 4502.64 16348 0 0 0 1 0 0 0 0 0 0 0
29 4356.98 15928 0 0 0 0 1 0 0 0 0 0 0
30 4591.27 16171 0 0 0 0 0 1 0 0 0 0 0
31 4696.96 15937 0 0 0 0 0 0 1 0 0 0 0
32 4621.40 15713 0 0 0 0 0 0 0 1 0 0 0
33 4562.84 15594 0 0 0 0 0 0 0 0 1 0 0
34 4202.52 15683 0 0 0 0 0 0 0 0 0 1 0
35 4296.49 16438 0 0 0 0 0 0 0 0 0 0 1
36 4435.23 17032 0 0 0 0 0 0 0 0 0 0 0
37 4105.18 17696 1 0 0 0 0 0 0 0 0 0 0
38 4116.68 17745 0 1 0 0 0 0 0 0 0 0 0
39 3844.49 19394 0 0 1 0 0 0 0 0 0 0 0
40 3720.98 20148 0 0 0 1 0 0 0 0 0 0 0
41 3674.40 20108 0 0 0 0 1 0 0 0 0 0 0
42 3857.62 18584 0 0 0 0 0 1 0 0 0 0 0
43 3801.06 18441 0 0 0 0 0 0 1 0 0 0 0
44 3504.37 18391 0 0 0 0 0 0 0 1 0 0 0
45 3032.60 19178 0 0 0 0 0 0 0 0 1 0 0
46 3047.03 18079 0 0 0 0 0 0 0 0 0 1 0
47 2962.34 18483 0 0 0 0 0 0 0 0 0 0 1
48 2197.82 19644 0 0 0 0 0 0 0 0 0 0 0
49 2014.45 19195 1 0 0 0 0 0 0 0 0 0 0
50 1862.83 19650 0 1 0 0 0 0 0 0 0 0 0
51 1905.41 20830 0 0 1 0 0 0 0 0 0 0 0
52 1810.99 23595 0 0 0 1 0 0 0 0 0 0 0
53 1670.07 22937 0 0 0 0 1 0 0 0 0 0 0
54 1864.44 21814 0 0 0 0 0 1 0 0 0 0 0
55 2052.02 21928 0 0 0 0 0 0 1 0 0 0 0
56 2029.60 21777 0 0 0 0 0 0 0 1 0 0 0
57 2070.83 21383 0 0 0 0 0 0 0 0 1 0 0
58 2293.41 21467 0 0 0 0 0 0 0 0 0 1 0
59 2443.27 22052 0 0 0 0 0 0 0 0 0 0 1
60 2513.17 22680 0 0 0 0 0 0 0 0 0 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
4956.68495 -0.09442 -206.66538 -162.86735 -63.22978 43.22000
M5 M6 M7 M8 M9 M10
-17.50933 80.42890 113.35963 -36.08779 -73.03666 -65.96042
M11
20.76740
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1113.05 -835.22 43.05 670.87 1184.38
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4956.68495 668.02912 7.420 1.89e-09 ***
X -0.09442 0.03121 -3.025 0.00402 **
M1 -206.66538 555.53854 -0.372 0.71156
M2 -162.86735 555.04636 -0.293 0.77049
M3 -63.22978 553.02740 -0.114 0.90946
M4 43.22000 551.77141 0.078 0.93790
M5 -17.50933 551.91750 -0.032 0.97483
M6 80.42890 552.17609 0.146 0.88481
M7 113.35963 551.98142 0.205 0.83817
M8 -36.08779 552.41109 -0.065 0.94819
M9 -73.03666 552.09439 -0.132 0.89532
M10 -65.96042 552.37459 -0.119 0.90546
M11 20.76740 551.98102 0.038 0.97015
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 872.3 on 47 degrees of freedom
Multiple R-squared: 0.1666, Adjusted R-squared: -0.04623
F-statistic: 0.7828 on 12 and 47 DF, p-value: 0.665
> 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,] 1.973955e-03 3.947911e-03 0.9980260
[2,] 2.501401e-04 5.002802e-04 0.9997499
[3,] 1.099021e-04 2.198043e-04 0.9998901
[4,] 2.125692e-04 4.251384e-04 0.9997874
[5,] 6.705591e-05 1.341118e-04 0.9999329
[6,] 2.327883e-05 4.655767e-05 0.9999767
[7,] 4.186218e-06 8.372436e-06 0.9999958
[8,] 3.295876e-06 6.591753e-06 0.9999967
[9,] 3.837499e-05 7.674998e-05 0.9999616
[10,] 1.883184e-04 3.766368e-04 0.9998117
[11,] 5.817697e-04 1.163539e-03 0.9994182
[12,] 1.128713e-03 2.257425e-03 0.9988713
[13,] 9.609301e-04 1.921860e-03 0.9990391
[14,] 6.590626e-04 1.318125e-03 0.9993409
[15,] 7.637531e-04 1.527506e-03 0.9992362
[16,] 2.713363e-03 5.426727e-03 0.9972866
[17,] 4.432263e-03 8.864526e-03 0.9955677
[18,] 5.226494e-03 1.045299e-02 0.9947735
[19,] 2.804499e-03 5.608997e-03 0.9971955
[20,] 1.283943e-03 2.567886e-03 0.9987161
[21,] 5.546020e-04 1.109204e-03 0.9994454
[22,] 1.188999e-03 2.377998e-03 0.9988110
[23,] 3.788762e-03 7.577523e-03 0.9962112
[24,] 3.662708e-02 7.325416e-02 0.9633729
[25,] 8.009253e-02 1.601851e-01 0.9199075
[26,] 1.812426e-01 3.624851e-01 0.8187574
[27,] 2.754251e-01 5.508502e-01 0.7245749
[28,] 3.573185e-01 7.146370e-01 0.6426815
[29,] 4.789913e-01 9.579825e-01 0.5210087
> postscript(file="/var/www/html/rcomp/tmp/1rctp1261305197.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/2qfg21261305197.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/3n13f1261305197.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/4tm711261305197.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/5xuq01261305197.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 = 60
Frequency = 1
1 2 3 4 5 6
-874.24318 -869.19231 -929.00288 -931.43529 -836.28513 -912.63580
7 8 9 10 11 12
-998.32915 -751.36597 -650.69802 -552.39043 -570.83211 -474.16007
13 14 15 16 17 18
-155.67593 4.44607 146.78206 222.90183 384.28238 369.54726
19 20 21 22 23 24
342.15459 81.64655 333.36161 467.02329 421.84442 593.75343
25 26 27 28 29 30
927.17451 942.13131 1021.38411 1046.27134 921.68535 1080.98055
31 32 33 34 35 36
1131.64614 1184.38406 1151.53726 792.54417 871.07151 1086.66287
37 38 39 40 41 42
1025.97143 998.29986 782.16664 623.39761 633.77024 575.15983
43 44 45 46 47 48
472.16741 320.20396 -40.31064 -136.72165 -269.99483 -904.12878
49 50 51 52 53 54
-923.22683 -1075.68492 -1021.32992 -961.13548 -1103.45283 -1113.05184
55 56 57 58 59 60
-947.63899 -834.86859 -793.89019 -570.45537 -452.08899 -302.12744
> postscript(file="/var/www/html/rcomp/tmp/6ogpt1261305197.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -874.24318 NA
1 -869.19231 -874.24318
2 -929.00288 -869.19231
3 -931.43529 -929.00288
4 -836.28513 -931.43529
5 -912.63580 -836.28513
6 -998.32915 -912.63580
7 -751.36597 -998.32915
8 -650.69802 -751.36597
9 -552.39043 -650.69802
10 -570.83211 -552.39043
11 -474.16007 -570.83211
12 -155.67593 -474.16007
13 4.44607 -155.67593
14 146.78206 4.44607
15 222.90183 146.78206
16 384.28238 222.90183
17 369.54726 384.28238
18 342.15459 369.54726
19 81.64655 342.15459
20 333.36161 81.64655
21 467.02329 333.36161
22 421.84442 467.02329
23 593.75343 421.84442
24 927.17451 593.75343
25 942.13131 927.17451
26 1021.38411 942.13131
27 1046.27134 1021.38411
28 921.68535 1046.27134
29 1080.98055 921.68535
30 1131.64614 1080.98055
31 1184.38406 1131.64614
32 1151.53726 1184.38406
33 792.54417 1151.53726
34 871.07151 792.54417
35 1086.66287 871.07151
36 1025.97143 1086.66287
37 998.29986 1025.97143
38 782.16664 998.29986
39 623.39761 782.16664
40 633.77024 623.39761
41 575.15983 633.77024
42 472.16741 575.15983
43 320.20396 472.16741
44 -40.31064 320.20396
45 -136.72165 -40.31064
46 -269.99483 -136.72165
47 -904.12878 -269.99483
48 -923.22683 -904.12878
49 -1075.68492 -923.22683
50 -1021.32992 -1075.68492
51 -961.13548 -1021.32992
52 -1103.45283 -961.13548
53 -1113.05184 -1103.45283
54 -947.63899 -1113.05184
55 -834.86859 -947.63899
56 -793.89019 -834.86859
57 -570.45537 -793.89019
58 -452.08899 -570.45537
59 -302.12744 -452.08899
60 NA -302.12744
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -869.19231 -874.24318
[2,] -929.00288 -869.19231
[3,] -931.43529 -929.00288
[4,] -836.28513 -931.43529
[5,] -912.63580 -836.28513
[6,] -998.32915 -912.63580
[7,] -751.36597 -998.32915
[8,] -650.69802 -751.36597
[9,] -552.39043 -650.69802
[10,] -570.83211 -552.39043
[11,] -474.16007 -570.83211
[12,] -155.67593 -474.16007
[13,] 4.44607 -155.67593
[14,] 146.78206 4.44607
[15,] 222.90183 146.78206
[16,] 384.28238 222.90183
[17,] 369.54726 384.28238
[18,] 342.15459 369.54726
[19,] 81.64655 342.15459
[20,] 333.36161 81.64655
[21,] 467.02329 333.36161
[22,] 421.84442 467.02329
[23,] 593.75343 421.84442
[24,] 927.17451 593.75343
[25,] 942.13131 927.17451
[26,] 1021.38411 942.13131
[27,] 1046.27134 1021.38411
[28,] 921.68535 1046.27134
[29,] 1080.98055 921.68535
[30,] 1131.64614 1080.98055
[31,] 1184.38406 1131.64614
[32,] 1151.53726 1184.38406
[33,] 792.54417 1151.53726
[34,] 871.07151 792.54417
[35,] 1086.66287 871.07151
[36,] 1025.97143 1086.66287
[37,] 998.29986 1025.97143
[38,] 782.16664 998.29986
[39,] 623.39761 782.16664
[40,] 633.77024 623.39761
[41,] 575.15983 633.77024
[42,] 472.16741 575.15983
[43,] 320.20396 472.16741
[44,] -40.31064 320.20396
[45,] -136.72165 -40.31064
[46,] -269.99483 -136.72165
[47,] -904.12878 -269.99483
[48,] -923.22683 -904.12878
[49,] -1075.68492 -923.22683
[50,] -1021.32992 -1075.68492
[51,] -961.13548 -1021.32992
[52,] -1103.45283 -961.13548
[53,] -1113.05184 -1103.45283
[54,] -947.63899 -1113.05184
[55,] -834.86859 -947.63899
[56,] -793.89019 -834.86859
[57,] -570.45537 -793.89019
[58,] -452.08899 -570.45537
[59,] -302.12744 -452.08899
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -869.19231 -874.24318
2 -929.00288 -869.19231
3 -931.43529 -929.00288
4 -836.28513 -931.43529
5 -912.63580 -836.28513
6 -998.32915 -912.63580
7 -751.36597 -998.32915
8 -650.69802 -751.36597
9 -552.39043 -650.69802
10 -570.83211 -552.39043
11 -474.16007 -570.83211
12 -155.67593 -474.16007
13 4.44607 -155.67593
14 146.78206 4.44607
15 222.90183 146.78206
16 384.28238 222.90183
17 369.54726 384.28238
18 342.15459 369.54726
19 81.64655 342.15459
20 333.36161 81.64655
21 467.02329 333.36161
22 421.84442 467.02329
23 593.75343 421.84442
24 927.17451 593.75343
25 942.13131 927.17451
26 1021.38411 942.13131
27 1046.27134 1021.38411
28 921.68535 1046.27134
29 1080.98055 921.68535
30 1131.64614 1080.98055
31 1184.38406 1131.64614
32 1151.53726 1184.38406
33 792.54417 1151.53726
34 871.07151 792.54417
35 1086.66287 871.07151
36 1025.97143 1086.66287
37 998.29986 1025.97143
38 782.16664 998.29986
39 623.39761 782.16664
40 633.77024 623.39761
41 575.15983 633.77024
42 472.16741 575.15983
43 320.20396 472.16741
44 -40.31064 320.20396
45 -136.72165 -40.31064
46 -269.99483 -136.72165
47 -904.12878 -269.99483
48 -923.22683 -904.12878
49 -1075.68492 -923.22683
50 -1021.32992 -1075.68492
51 -961.13548 -1021.32992
52 -1103.45283 -961.13548
53 -1113.05184 -1103.45283
54 -947.63899 -1113.05184
55 -834.86859 -947.63899
56 -793.89019 -834.86859
57 -570.45537 -793.89019
58 -452.08899 -570.45537
59 -302.12744 -452.08899
> 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/7c8y01261305197.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/8gd6v1261305197.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/91w7t1261305197.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/10c31i1261305197.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/11djg41261305197.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/12gjn81261305197.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/13lzt41261305197.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/1493041261305197.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/15h30c1261305197.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/16d5kf1261305197.tab")
+ }
>
> try(system("convert tmp/1rctp1261305197.ps tmp/1rctp1261305197.png",intern=TRUE))
character(0)
> try(system("convert tmp/2qfg21261305197.ps tmp/2qfg21261305197.png",intern=TRUE))
character(0)
> try(system("convert tmp/3n13f1261305197.ps tmp/3n13f1261305197.png",intern=TRUE))
character(0)
> try(system("convert tmp/4tm711261305197.ps tmp/4tm711261305197.png",intern=TRUE))
character(0)
> try(system("convert tmp/5xuq01261305197.ps tmp/5xuq01261305197.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ogpt1261305197.ps tmp/6ogpt1261305197.png",intern=TRUE))
character(0)
> try(system("convert tmp/7c8y01261305197.ps tmp/7c8y01261305197.png",intern=TRUE))
character(0)
> try(system("convert tmp/8gd6v1261305197.ps tmp/8gd6v1261305197.png",intern=TRUE))
character(0)
> try(system("convert tmp/91w7t1261305197.ps tmp/91w7t1261305197.png",intern=TRUE))
character(0)
> try(system("convert tmp/10c31i1261305197.ps tmp/10c31i1261305197.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.375 1.562 3.560