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 = '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 2849.27 10872 1 0 0 0 0 0 0 0 0 0 0 1
2 2921.44 10625 0 1 0 0 0 0 0 0 0 0 0 2
3 2981.85 10407 0 0 1 0 0 0 0 0 0 0 0 3
4 3080.58 10463 0 0 0 1 0 0 0 0 0 0 0 4
5 3106.22 10556 0 0 0 0 1 0 0 0 0 0 0 5
6 3119.31 10646 0 0 0 0 0 1 0 0 0 0 0 6
7 3061.26 10702 0 0 0 0 0 0 1 0 0 0 0 7
8 3097.31 11353 0 0 0 0 0 0 0 1 0 0 0 8
9 3161.69 11346 0 0 0 0 0 0 0 0 1 0 0 9
10 3257.16 11451 0 0 0 0 0 0 0 0 0 1 0 10
11 3277.01 11964 0 0 0 0 0 0 0 0 0 0 1 11
12 3295.32 12574 0 0 0 0 0 0 0 0 0 0 0 12
13 3363.99 13031 1 0 0 0 0 0 0 0 0 0 0 13
14 3494.17 13812 0 1 0 0 0 0 0 0 0 0 0 14
15 3667.03 14544 0 0 1 0 0 0 0 0 0 0 0 15
16 3813.06 14931 0 0 0 1 0 0 0 0 0 0 0 16
17 3917.96 14886 0 0 0 0 1 0 0 0 0 0 0 17
18 3895.51 16005 0 0 0 0 0 1 0 0 0 0 0 18
19 3801.06 17064 0 0 0 0 0 0 1 0 0 0 0 19
20 3570.12 15168 0 0 0 0 0 0 0 1 0 0 0 20
21 3701.61 16050 0 0 0 0 0 0 0 0 1 0 0 21
22 3862.27 15839 0 0 0 0 0 0 0 0 0 1 0 22
23 3970.10 15137 0 0 0 0 0 0 0 0 0 0 1 23
24 4138.52 14954 0 0 0 0 0 0 0 0 0 0 0 24
25 4199.75 15648 1 0 0 0 0 0 0 0 0 0 0 25
26 4290.89 15305 0 1 0 0 0 0 0 0 0 0 0 26
27 4443.91 15579 0 0 1 0 0 0 0 0 0 0 0 27
28 4502.64 16348 0 0 0 1 0 0 0 0 0 0 0 28
29 4356.98 15928 0 0 0 0 1 0 0 0 0 0 0 29
30 4591.27 16171 0 0 0 0 0 1 0 0 0 0 0 30
31 4696.96 15937 0 0 0 0 0 0 1 0 0 0 0 31
32 4621.40 15713 0 0 0 0 0 0 0 1 0 0 0 32
33 4562.84 15594 0 0 0 0 0 0 0 0 1 0 0 33
34 4202.52 15683 0 0 0 0 0 0 0 0 0 1 0 34
35 4296.49 16438 0 0 0 0 0 0 0 0 0 0 1 35
36 4435.23 17032 0 0 0 0 0 0 0 0 0 0 0 36
37 4105.18 17696 1 0 0 0 0 0 0 0 0 0 0 37
38 4116.68 17745 0 1 0 0 0 0 0 0 0 0 0 38
39 3844.49 19394 0 0 1 0 0 0 0 0 0 0 0 39
40 3720.98 20148 0 0 0 1 0 0 0 0 0 0 0 40
41 3674.40 20108 0 0 0 0 1 0 0 0 0 0 0 41
42 3857.62 18584 0 0 0 0 0 1 0 0 0 0 0 42
43 3801.06 18441 0 0 0 0 0 0 1 0 0 0 0 43
44 3504.37 18391 0 0 0 0 0 0 0 1 0 0 0 44
45 3032.60 19178 0 0 0 0 0 0 0 0 1 0 0 45
46 3047.03 18079 0 0 0 0 0 0 0 0 0 1 0 46
47 2962.34 18483 0 0 0 0 0 0 0 0 0 0 1 47
48 2197.82 19644 0 0 0 0 0 0 0 0 0 0 0 48
49 2014.45 19195 1 0 0 0 0 0 0 0 0 0 0 49
50 1862.83 19650 0 1 0 0 0 0 0 0 0 0 0 50
51 1905.41 20830 0 0 1 0 0 0 0 0 0 0 0 51
52 1810.99 23595 0 0 0 1 0 0 0 0 0 0 0 52
53 1670.07 22937 0 0 0 0 1 0 0 0 0 0 0 53
54 1864.44 21814 0 0 0 0 0 1 0 0 0 0 0 54
55 2052.02 21928 0 0 0 0 0 0 1 0 0 0 0 55
56 2029.60 21777 0 0 0 0 0 0 0 1 0 0 0 56
57 2070.83 21383 0 0 0 0 0 0 0 0 1 0 0 57
58 2293.41 21467 0 0 0 0 0 0 0 0 0 1 0 58
59 2443.27 22052 0 0 0 0 0 0 0 0 0 0 1 59
60 2513.17 22680 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
4631.31499 -0.06013 -217.69707 -171.15242 -88.80660 -7.28811
M5 M6 M7 M8 M9 M10
-53.16687 60.47914 95.07983 -35.40238 -72.71805 -51.05184
M11 t
32.52494 -7.51276
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1130.731 -808.779 3.089 650.285 1210.697
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4631.31499 1380.13614 3.356 0.00159 **
X -0.06013 0.13074 -0.460 0.64774
M1 -217.69707 562.58220 -0.387 0.70057
M2 -171.15242 561.43982 -0.305 0.76186
M3 -88.80660 566.52425 -0.157 0.87612
M4 -7.28811 587.79833 -0.012 0.99016
M5 -53.16687 572.84438 -0.093 0.92646
M6 60.47914 562.56754 0.108 0.91486
M7 95.07983 561.59490 0.169 0.86630
M8 -35.40238 557.94630 -0.063 0.94968
M9 -72.71805 557.62190 -0.130 0.89681
M10 -51.05184 560.62449 -0.091 0.92784
M11 32.52494 559.20115 0.058 0.95387
t -7.51276 27.79947 -0.270 0.78818
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 881 on 46 degrees of freedom
Multiple R-squared: 0.1679, Adjusted R-squared: -0.06727
F-statistic: 0.7139 on 13 and 46 DF, p-value: 0.7401
> 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.887685e-04 9.775370e-04 0.9995112
[2,] 5.322138e-05 1.064428e-04 0.9999468
[3,] 2.192480e-05 4.384960e-05 0.9999781
[4,] 2.321764e-05 4.643528e-05 0.9999768
[5,] 9.072652e-06 1.814530e-05 0.9999909
[6,] 1.533703e-06 3.067405e-06 0.9999985
[7,] 7.783089e-07 1.556618e-06 0.9999992
[8,] 3.913185e-06 7.826370e-06 0.9999961
[9,] 2.236797e-06 4.473594e-06 0.9999978
[10,] 9.049633e-07 1.809927e-06 0.9999991
[11,] 2.368077e-07 4.736153e-07 0.9999998
[12,] 4.359014e-08 8.718028e-08 1.0000000
[13,] 3.955055e-08 7.910110e-08 1.0000000
[14,] 9.075595e-09 1.815119e-08 1.0000000
[15,] 9.077113e-09 1.815423e-08 1.0000000
[16,] 6.160524e-09 1.232105e-08 1.0000000
[17,] 1.205285e-09 2.410571e-09 1.0000000
[18,] 4.955040e-08 9.910080e-08 1.0000000
[19,] 2.146359e-07 4.292717e-07 0.9999998
[20,] 1.757909e-07 3.515817e-07 0.9999998
[21,] 5.908267e-06 1.181653e-05 0.9999941
[22,] 6.024603e-05 1.204921e-04 0.9999398
[23,] 1.257898e-03 2.515796e-03 0.9987421
[24,] 8.627996e-03 1.725599e-02 0.9913720
[25,] 3.288578e-02 6.577156e-02 0.9671142
[26,] 9.233515e-02 1.846703e-01 0.9076648
[27,] 1.757113e-01 3.514225e-01 0.8242887
> postscript(file="/var/www/html/rcomp/tmp/1jiuu1261313776.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/2q4v21261313776.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/3gep91261313776.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/43ack1261313776.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/5h0d91261313776.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
-903.11649 -884.83015 -912.36126 -884.26978 -799.64629 -887.27796
7 8 9 10 11 12
-969.04868 -755.85996 -647.07243 -559.44236 -584.81038 -489.78420
13 14 15 16 17 18
-168.42558 -30.31698 111.72412 207.01818 362.60392 301.30463
19 20 21 22 23 24
243.44295 36.49400 265.84590 399.66531 389.22099 586.67515
25 26 27 28 29 30
914.84426 946.32825 1040.99044 1071.95364 954.43114 1097.19915
31 32 33 34 35 36
1161.73112 1210.69728 1189.81040 820.68841 883.99153 1098.48564
37 38 39 40 41 42
1033.57090 1008.98532 761.11440 608.93567 613.34206 598.79274
43 44 45 46 47 48
506.54642 344.84496 -34.77536 -100.58019 -237.04222 -891.71516
49 50 51 52 53 54
-876.87309 -1040.16644 -1001.46770 -1003.63771 -1130.73083 -1110.01856
55 56 57 58 59 60
-942.67181 -836.17627 -773.80852 -560.33116 -451.35991 -303.66142
> postscript(file="/var/www/html/rcomp/tmp/6zo7h1261313776.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 -903.11649 NA
1 -884.83015 -903.11649
2 -912.36126 -884.83015
3 -884.26978 -912.36126
4 -799.64629 -884.26978
5 -887.27796 -799.64629
6 -969.04868 -887.27796
7 -755.85996 -969.04868
8 -647.07243 -755.85996
9 -559.44236 -647.07243
10 -584.81038 -559.44236
11 -489.78420 -584.81038
12 -168.42558 -489.78420
13 -30.31698 -168.42558
14 111.72412 -30.31698
15 207.01818 111.72412
16 362.60392 207.01818
17 301.30463 362.60392
18 243.44295 301.30463
19 36.49400 243.44295
20 265.84590 36.49400
21 399.66531 265.84590
22 389.22099 399.66531
23 586.67515 389.22099
24 914.84426 586.67515
25 946.32825 914.84426
26 1040.99044 946.32825
27 1071.95364 1040.99044
28 954.43114 1071.95364
29 1097.19915 954.43114
30 1161.73112 1097.19915
31 1210.69728 1161.73112
32 1189.81040 1210.69728
33 820.68841 1189.81040
34 883.99153 820.68841
35 1098.48564 883.99153
36 1033.57090 1098.48564
37 1008.98532 1033.57090
38 761.11440 1008.98532
39 608.93567 761.11440
40 613.34206 608.93567
41 598.79274 613.34206
42 506.54642 598.79274
43 344.84496 506.54642
44 -34.77536 344.84496
45 -100.58019 -34.77536
46 -237.04222 -100.58019
47 -891.71516 -237.04222
48 -876.87309 -891.71516
49 -1040.16644 -876.87309
50 -1001.46770 -1040.16644
51 -1003.63771 -1001.46770
52 -1130.73083 -1003.63771
53 -1110.01856 -1130.73083
54 -942.67181 -1110.01856
55 -836.17627 -942.67181
56 -773.80852 -836.17627
57 -560.33116 -773.80852
58 -451.35991 -560.33116
59 -303.66142 -451.35991
60 NA -303.66142
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -884.83015 -903.11649
[2,] -912.36126 -884.83015
[3,] -884.26978 -912.36126
[4,] -799.64629 -884.26978
[5,] -887.27796 -799.64629
[6,] -969.04868 -887.27796
[7,] -755.85996 -969.04868
[8,] -647.07243 -755.85996
[9,] -559.44236 -647.07243
[10,] -584.81038 -559.44236
[11,] -489.78420 -584.81038
[12,] -168.42558 -489.78420
[13,] -30.31698 -168.42558
[14,] 111.72412 -30.31698
[15,] 207.01818 111.72412
[16,] 362.60392 207.01818
[17,] 301.30463 362.60392
[18,] 243.44295 301.30463
[19,] 36.49400 243.44295
[20,] 265.84590 36.49400
[21,] 399.66531 265.84590
[22,] 389.22099 399.66531
[23,] 586.67515 389.22099
[24,] 914.84426 586.67515
[25,] 946.32825 914.84426
[26,] 1040.99044 946.32825
[27,] 1071.95364 1040.99044
[28,] 954.43114 1071.95364
[29,] 1097.19915 954.43114
[30,] 1161.73112 1097.19915
[31,] 1210.69728 1161.73112
[32,] 1189.81040 1210.69728
[33,] 820.68841 1189.81040
[34,] 883.99153 820.68841
[35,] 1098.48564 883.99153
[36,] 1033.57090 1098.48564
[37,] 1008.98532 1033.57090
[38,] 761.11440 1008.98532
[39,] 608.93567 761.11440
[40,] 613.34206 608.93567
[41,] 598.79274 613.34206
[42,] 506.54642 598.79274
[43,] 344.84496 506.54642
[44,] -34.77536 344.84496
[45,] -100.58019 -34.77536
[46,] -237.04222 -100.58019
[47,] -891.71516 -237.04222
[48,] -876.87309 -891.71516
[49,] -1040.16644 -876.87309
[50,] -1001.46770 -1040.16644
[51,] -1003.63771 -1001.46770
[52,] -1130.73083 -1003.63771
[53,] -1110.01856 -1130.73083
[54,] -942.67181 -1110.01856
[55,] -836.17627 -942.67181
[56,] -773.80852 -836.17627
[57,] -560.33116 -773.80852
[58,] -451.35991 -560.33116
[59,] -303.66142 -451.35991
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -884.83015 -903.11649
2 -912.36126 -884.83015
3 -884.26978 -912.36126
4 -799.64629 -884.26978
5 -887.27796 -799.64629
6 -969.04868 -887.27796
7 -755.85996 -969.04868
8 -647.07243 -755.85996
9 -559.44236 -647.07243
10 -584.81038 -559.44236
11 -489.78420 -584.81038
12 -168.42558 -489.78420
13 -30.31698 -168.42558
14 111.72412 -30.31698
15 207.01818 111.72412
16 362.60392 207.01818
17 301.30463 362.60392
18 243.44295 301.30463
19 36.49400 243.44295
20 265.84590 36.49400
21 399.66531 265.84590
22 389.22099 399.66531
23 586.67515 389.22099
24 914.84426 586.67515
25 946.32825 914.84426
26 1040.99044 946.32825
27 1071.95364 1040.99044
28 954.43114 1071.95364
29 1097.19915 954.43114
30 1161.73112 1097.19915
31 1210.69728 1161.73112
32 1189.81040 1210.69728
33 820.68841 1189.81040
34 883.99153 820.68841
35 1098.48564 883.99153
36 1033.57090 1098.48564
37 1008.98532 1033.57090
38 761.11440 1008.98532
39 608.93567 761.11440
40 613.34206 608.93567
41 598.79274 613.34206
42 506.54642 598.79274
43 344.84496 506.54642
44 -34.77536 344.84496
45 -100.58019 -34.77536
46 -237.04222 -100.58019
47 -891.71516 -237.04222
48 -876.87309 -891.71516
49 -1040.16644 -876.87309
50 -1001.46770 -1040.16644
51 -1003.63771 -1001.46770
52 -1130.73083 -1003.63771
53 -1110.01856 -1130.73083
54 -942.67181 -1110.01856
55 -836.17627 -942.67181
56 -773.80852 -836.17627
57 -560.33116 -773.80852
58 -451.35991 -560.33116
59 -303.66142 -451.35991
> 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/7nb6v1261313776.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/855161261313776.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/9wc5n1261313776.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/10o8db1261313776.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/1138jj1261313776.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/12pu2k1261313776.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/136nng1261313776.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/144n9w1261313776.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/152kdp1261313776.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/16rf2h1261313776.tab")
+ }
>
> try(system("convert tmp/1jiuu1261313776.ps tmp/1jiuu1261313776.png",intern=TRUE))
character(0)
> try(system("convert tmp/2q4v21261313776.ps tmp/2q4v21261313776.png",intern=TRUE))
character(0)
> try(system("convert tmp/3gep91261313776.ps tmp/3gep91261313776.png",intern=TRUE))
character(0)
> try(system("convert tmp/43ack1261313776.ps tmp/43ack1261313776.png",intern=TRUE))
character(0)
> try(system("convert tmp/5h0d91261313776.ps tmp/5h0d91261313776.png",intern=TRUE))
character(0)
> try(system("convert tmp/6zo7h1261313776.ps tmp/6zo7h1261313776.png",intern=TRUE))
character(0)
> try(system("convert tmp/7nb6v1261313776.ps tmp/7nb6v1261313776.png",intern=TRUE))
character(0)
> try(system("convert tmp/855161261313776.ps tmp/855161261313776.png",intern=TRUE))
character(0)
> try(system("convert tmp/9wc5n1261313776.ps tmp/9wc5n1261313776.png",intern=TRUE))
character(0)
> try(system("convert tmp/10o8db1261313776.ps tmp/10o8db1261313776.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.325 1.583 3.324