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(269285,8.2,269829,8,270911,7.5,266844,6.8,271244,6.5,269907,6.6,271296,7.6,270157,8,271322,8.1,267179,7.7,264101,7.5,265518,7.6,269419,7.8,268714,7.8,272482,7.8,268351,7.5,268175,7.5,270674,7.1,272764,7.5,272599,7.5,270333,7.6,270846,7.7,270491,7.7,269160,7.9,274027,8.1,273784,8.2,276663,8.2,274525,8.2,271344,7.9,271115,7.3,270798,6.9,273911,6.6,273985,6.7,271917,6.9,273338,7,270601,7.1,273547,7.2,275363,7.1,281229,6.9,277793,7,279913,6.8,282500,6.4,280041,6.7,282166,6.6,290304,6.4,283519,6.3,287816,6.2,285226,6.5,287595,6.8,289741,6.8,289148,6.4,288301,6.1,290155,5.8,289648,6.1,288225,7.2,289351,7.3,294735,6.9,305333,6.1,309030,5.8,310215,6.2,321935,7.1,325734,7.7,320846,7.9,323023,7.7,319753,7.4,321753,7.5,320757,8,324479,8.1),dim=c(2,68),dimnames=list(c('Y','X'),1:68))
> y <- array(NA,dim=c(2,68),dimnames=list(c('Y','X'),1:68))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> #'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
1 269285 8.2
2 269829 8.0
3 270911 7.5
4 266844 6.8
5 271244 6.5
6 269907 6.6
7 271296 7.6
8 270157 8.0
9 271322 8.1
10 267179 7.7
11 264101 7.5
12 265518 7.6
13 269419 7.8
14 268714 7.8
15 272482 7.8
16 268351 7.5
17 268175 7.5
18 270674 7.1
19 272764 7.5
20 272599 7.5
21 270333 7.6
22 270846 7.7
23 270491 7.7
24 269160 7.9
25 274027 8.1
26 273784 8.2
27 276663 8.2
28 274525 8.2
29 271344 7.9
30 271115 7.3
31 270798 6.9
32 273911 6.6
33 273985 6.7
34 271917 6.9
35 273338 7.0
36 270601 7.1
37 273547 7.2
38 275363 7.1
39 281229 6.9
40 277793 7.0
41 279913 6.8
42 282500 6.4
43 280041 6.7
44 282166 6.6
45 290304 6.4
46 283519 6.3
47 287816 6.2
48 285226 6.5
49 287595 6.8
50 289741 6.8
51 289148 6.4
52 288301 6.1
53 290155 5.8
54 289648 6.1
55 288225 7.2
56 289351 7.3
57 294735 6.9
58 305333 6.1
59 309030 5.8
60 310215 6.2
61 321935 7.1
62 325734 7.7
63 320846 7.9
64 323023 7.7
65 319753 7.4
66 321753 7.5
67 320757 8.0
68 324479 8.1
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X
311563 -3982
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-17643 -10952 -6913 3079 45169
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 311563 23352 13.342 <2e-16 ***
X -3982 3230 -1.233 0.222
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 17480 on 66 degrees of freedom
Multiple R-squared: 0.02251, Adjusted R-squared: 0.007697
F-statistic: 1.52 on 1 and 66 DF, p-value: 0.2220
> 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,] 2.005849e-03 4.011697e-03 0.9979941513
[2,] 1.901582e-04 3.803163e-04 0.9998098418
[3,] 2.401511e-05 4.803022e-05 0.9999759849
[4,] 2.015992e-06 4.031983e-06 0.9999979840
[5,] 2.077170e-07 4.154339e-07 0.9999997923
[6,] 8.366756e-08 1.673351e-07 0.9999999163
[7,] 2.904165e-07 5.808329e-07 0.9999997096
[8,] 1.160862e-07 2.321724e-07 0.9999998839
[9,] 1.651569e-08 3.303137e-08 0.9999999835
[10,] 2.292592e-09 4.585184e-09 0.9999999977
[11,] 7.854984e-10 1.570997e-09 0.9999999992
[12,] 1.192136e-10 2.384272e-10 0.9999999999
[13,] 1.842680e-11 3.685361e-11 1.0000000000
[14,] 3.284804e-12 6.569609e-12 1.0000000000
[15,] 1.372756e-12 2.745511e-12 1.0000000000
[16,] 4.627426e-13 9.254852e-13 1.0000000000
[17,] 7.446659e-14 1.489332e-13 1.0000000000
[18,] 1.304890e-14 2.609781e-14 1.0000000000
[19,] 2.197935e-15 4.395869e-15 1.0000000000
[20,] 3.976417e-16 7.952834e-16 1.0000000000
[21,] 3.219396e-16 6.438792e-16 1.0000000000
[22,] 1.730589e-16 3.461179e-16 1.0000000000
[23,] 5.058579e-16 1.011716e-15 1.0000000000
[24,] 3.322441e-16 6.644883e-16 1.0000000000
[25,] 1.606757e-16 3.213514e-16 1.0000000000
[26,] 7.140702e-17 1.428140e-16 1.0000000000
[27,] 2.850257e-17 5.700514e-17 1.0000000000
[28,] 3.633510e-17 7.267019e-17 1.0000000000
[29,] 3.107330e-17 6.214660e-17 1.0000000000
[30,] 1.566495e-17 3.132990e-17 1.0000000000
[31,] 1.344046e-17 2.688092e-17 1.0000000000
[32,] 1.432471e-17 2.864942e-17 1.0000000000
[33,] 3.418502e-17 6.837004e-17 1.0000000000
[34,] 1.847433e-16 3.694866e-16 1.0000000000
[35,] 2.286382e-14 4.572763e-14 1.0000000000
[36,] 2.120385e-13 4.240771e-13 1.0000000000
[37,] 2.318539e-12 4.637078e-12 1.0000000000
[38,] 1.426167e-11 2.852333e-11 1.0000000000
[39,] 6.679737e-11 1.335947e-10 0.9999999999
[40,] 3.435038e-10 6.870076e-10 0.9999999997
[41,] 9.544790e-09 1.908958e-08 0.9999999905
[42,] 1.204322e-08 2.408645e-08 0.9999999880
[43,] 2.027522e-08 4.055044e-08 0.9999999797
[44,] 4.582740e-08 9.165480e-08 0.9999999542
[45,] 3.874873e-07 7.749747e-07 0.9999996125
[46,] 3.460524e-06 6.921049e-06 0.9999965395
[47,] 6.923245e-06 1.384649e-05 0.9999930768
[48,] 7.112126e-06 1.422425e-05 0.9999928879
[49,] 4.631223e-06 9.262445e-06 0.9999953688
[50,] 6.824537e-06 1.364907e-05 0.9999931755
[51,] 8.505549e-04 1.701110e-03 0.9991494451
[52,] 1.953930e-01 3.907859e-01 0.8046070257
[53,] 9.888490e-01 2.230192e-02 0.0111509589
[54,] 9.966572e-01 6.685638e-03 0.0033428191
[55,] 9.941179e-01 1.176423e-02 0.0058821154
[56,] 9.982219e-01 3.556150e-03 0.0017780748
[57,] 9.977329e-01 4.534173e-03 0.0022670863
[58,] 9.991382e-01 1.723599e-03 0.0008617996
[59,] 9.968337e-01 6.332685e-03 0.0031663425
> postscript(file="/var/www/html/rcomp/tmp/1fj2n1258804002.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/2heut1258804002.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/3cylv1258804002.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/4etd01258804002.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/58fke1258804002.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 = 68
Frequency = 1
1 2 3 4 5 6
-9627.0917 -9879.4469 -10788.3349 -17642.5780 -14437.1108 -15375.9332
7 8 9 10 11 12
-10005.1573 -9551.4469 -7988.2693 -13723.9797 -17598.3349 -15783.1573
13 14 15 16 17 18
-11085.8021 -11790.8021 -8022.8021 -13348.3349 -13524.3349 -12618.0452
19 20 21 22 23 24
-8935.3349 -9100.3349 -10968.1573 -10056.9797 -10411.9797 -10946.6245
25 26 27 28 29 30
-5283.2693 -5128.0917 -2249.0917 -4387.0917 -8762.6245 -11380.6901
31 32 33 34 35 36
-13290.4004 -11371.9332 -10899.7556 -12171.4004 -10352.2228 -12691.0452
37 38 39 40 41 42
-9346.8676 -7929.0452 -2859.4004 -5897.2228 -4573.5780 -3579.2884
43 44 45 46 47 48
-4843.7556 -3116.9332 4224.7116 -2958.4660 940.3564 -455.1108
49 50 51 52 53 54
3108.4220 5254.4220 3068.7116 1027.1788 1686.6460 2374.1788
55 56 57 58 59 60
5331.1324 6855.3099 10646.5996 18059.1788 20561.6460 23339.3564
61 62 63 64 65 66
38642.9548 44831.0203 40739.3755 42120.0203 37655.4875 40053.6651
67 68
41048.5531 45168.7307
> postscript(file="/var/www/html/rcomp/tmp/6jc4f1258804002.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 = 68
Frequency = 1
lag(myerror, k = 1) myerror
0 -9627.0917 NA
1 -9879.4469 -9627.0917
2 -10788.3349 -9879.4469
3 -17642.5780 -10788.3349
4 -14437.1108 -17642.5780
5 -15375.9332 -14437.1108
6 -10005.1573 -15375.9332
7 -9551.4469 -10005.1573
8 -7988.2693 -9551.4469
9 -13723.9797 -7988.2693
10 -17598.3349 -13723.9797
11 -15783.1573 -17598.3349
12 -11085.8021 -15783.1573
13 -11790.8021 -11085.8021
14 -8022.8021 -11790.8021
15 -13348.3349 -8022.8021
16 -13524.3349 -13348.3349
17 -12618.0452 -13524.3349
18 -8935.3349 -12618.0452
19 -9100.3349 -8935.3349
20 -10968.1573 -9100.3349
21 -10056.9797 -10968.1573
22 -10411.9797 -10056.9797
23 -10946.6245 -10411.9797
24 -5283.2693 -10946.6245
25 -5128.0917 -5283.2693
26 -2249.0917 -5128.0917
27 -4387.0917 -2249.0917
28 -8762.6245 -4387.0917
29 -11380.6901 -8762.6245
30 -13290.4004 -11380.6901
31 -11371.9332 -13290.4004
32 -10899.7556 -11371.9332
33 -12171.4004 -10899.7556
34 -10352.2228 -12171.4004
35 -12691.0452 -10352.2228
36 -9346.8676 -12691.0452
37 -7929.0452 -9346.8676
38 -2859.4004 -7929.0452
39 -5897.2228 -2859.4004
40 -4573.5780 -5897.2228
41 -3579.2884 -4573.5780
42 -4843.7556 -3579.2884
43 -3116.9332 -4843.7556
44 4224.7116 -3116.9332
45 -2958.4660 4224.7116
46 940.3564 -2958.4660
47 -455.1108 940.3564
48 3108.4220 -455.1108
49 5254.4220 3108.4220
50 3068.7116 5254.4220
51 1027.1788 3068.7116
52 1686.6460 1027.1788
53 2374.1788 1686.6460
54 5331.1324 2374.1788
55 6855.3099 5331.1324
56 10646.5996 6855.3099
57 18059.1788 10646.5996
58 20561.6460 18059.1788
59 23339.3564 20561.6460
60 38642.9548 23339.3564
61 44831.0203 38642.9548
62 40739.3755 44831.0203
63 42120.0203 40739.3755
64 37655.4875 42120.0203
65 40053.6651 37655.4875
66 41048.5531 40053.6651
67 45168.7307 41048.5531
68 NA 45168.7307
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -9879.4469 -9627.0917
[2,] -10788.3349 -9879.4469
[3,] -17642.5780 -10788.3349
[4,] -14437.1108 -17642.5780
[5,] -15375.9332 -14437.1108
[6,] -10005.1573 -15375.9332
[7,] -9551.4469 -10005.1573
[8,] -7988.2693 -9551.4469
[9,] -13723.9797 -7988.2693
[10,] -17598.3349 -13723.9797
[11,] -15783.1573 -17598.3349
[12,] -11085.8021 -15783.1573
[13,] -11790.8021 -11085.8021
[14,] -8022.8021 -11790.8021
[15,] -13348.3349 -8022.8021
[16,] -13524.3349 -13348.3349
[17,] -12618.0452 -13524.3349
[18,] -8935.3349 -12618.0452
[19,] -9100.3349 -8935.3349
[20,] -10968.1573 -9100.3349
[21,] -10056.9797 -10968.1573
[22,] -10411.9797 -10056.9797
[23,] -10946.6245 -10411.9797
[24,] -5283.2693 -10946.6245
[25,] -5128.0917 -5283.2693
[26,] -2249.0917 -5128.0917
[27,] -4387.0917 -2249.0917
[28,] -8762.6245 -4387.0917
[29,] -11380.6901 -8762.6245
[30,] -13290.4004 -11380.6901
[31,] -11371.9332 -13290.4004
[32,] -10899.7556 -11371.9332
[33,] -12171.4004 -10899.7556
[34,] -10352.2228 -12171.4004
[35,] -12691.0452 -10352.2228
[36,] -9346.8676 -12691.0452
[37,] -7929.0452 -9346.8676
[38,] -2859.4004 -7929.0452
[39,] -5897.2228 -2859.4004
[40,] -4573.5780 -5897.2228
[41,] -3579.2884 -4573.5780
[42,] -4843.7556 -3579.2884
[43,] -3116.9332 -4843.7556
[44,] 4224.7116 -3116.9332
[45,] -2958.4660 4224.7116
[46,] 940.3564 -2958.4660
[47,] -455.1108 940.3564
[48,] 3108.4220 -455.1108
[49,] 5254.4220 3108.4220
[50,] 3068.7116 5254.4220
[51,] 1027.1788 3068.7116
[52,] 1686.6460 1027.1788
[53,] 2374.1788 1686.6460
[54,] 5331.1324 2374.1788
[55,] 6855.3099 5331.1324
[56,] 10646.5996 6855.3099
[57,] 18059.1788 10646.5996
[58,] 20561.6460 18059.1788
[59,] 23339.3564 20561.6460
[60,] 38642.9548 23339.3564
[61,] 44831.0203 38642.9548
[62,] 40739.3755 44831.0203
[63,] 42120.0203 40739.3755
[64,] 37655.4875 42120.0203
[65,] 40053.6651 37655.4875
[66,] 41048.5531 40053.6651
[67,] 45168.7307 41048.5531
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -9879.4469 -9627.0917
2 -10788.3349 -9879.4469
3 -17642.5780 -10788.3349
4 -14437.1108 -17642.5780
5 -15375.9332 -14437.1108
6 -10005.1573 -15375.9332
7 -9551.4469 -10005.1573
8 -7988.2693 -9551.4469
9 -13723.9797 -7988.2693
10 -17598.3349 -13723.9797
11 -15783.1573 -17598.3349
12 -11085.8021 -15783.1573
13 -11790.8021 -11085.8021
14 -8022.8021 -11790.8021
15 -13348.3349 -8022.8021
16 -13524.3349 -13348.3349
17 -12618.0452 -13524.3349
18 -8935.3349 -12618.0452
19 -9100.3349 -8935.3349
20 -10968.1573 -9100.3349
21 -10056.9797 -10968.1573
22 -10411.9797 -10056.9797
23 -10946.6245 -10411.9797
24 -5283.2693 -10946.6245
25 -5128.0917 -5283.2693
26 -2249.0917 -5128.0917
27 -4387.0917 -2249.0917
28 -8762.6245 -4387.0917
29 -11380.6901 -8762.6245
30 -13290.4004 -11380.6901
31 -11371.9332 -13290.4004
32 -10899.7556 -11371.9332
33 -12171.4004 -10899.7556
34 -10352.2228 -12171.4004
35 -12691.0452 -10352.2228
36 -9346.8676 -12691.0452
37 -7929.0452 -9346.8676
38 -2859.4004 -7929.0452
39 -5897.2228 -2859.4004
40 -4573.5780 -5897.2228
41 -3579.2884 -4573.5780
42 -4843.7556 -3579.2884
43 -3116.9332 -4843.7556
44 4224.7116 -3116.9332
45 -2958.4660 4224.7116
46 940.3564 -2958.4660
47 -455.1108 940.3564
48 3108.4220 -455.1108
49 5254.4220 3108.4220
50 3068.7116 5254.4220
51 1027.1788 3068.7116
52 1686.6460 1027.1788
53 2374.1788 1686.6460
54 5331.1324 2374.1788
55 6855.3099 5331.1324
56 10646.5996 6855.3099
57 18059.1788 10646.5996
58 20561.6460 18059.1788
59 23339.3564 20561.6460
60 38642.9548 23339.3564
61 44831.0203 38642.9548
62 40739.3755 44831.0203
63 42120.0203 40739.3755
64 37655.4875 42120.0203
65 40053.6651 37655.4875
66 41048.5531 40053.6651
67 45168.7307 41048.5531
> 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/7gut21258804002.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/8zuez1258804002.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/99j691258804002.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/10nzyh1258804002.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/111hx41258804002.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/12b2bh1258804002.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/13qqgm1258804002.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/14xgcl1258804002.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/15p6ry1258804002.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/16gmzo1258804002.tab")
+ }
>
> system("convert tmp/1fj2n1258804002.ps tmp/1fj2n1258804002.png")
> system("convert tmp/2heut1258804002.ps tmp/2heut1258804002.png")
> system("convert tmp/3cylv1258804002.ps tmp/3cylv1258804002.png")
> system("convert tmp/4etd01258804002.ps tmp/4etd01258804002.png")
> system("convert tmp/58fke1258804002.ps tmp/58fke1258804002.png")
> system("convert tmp/6jc4f1258804002.ps tmp/6jc4f1258804002.png")
> system("convert tmp/7gut21258804002.ps tmp/7gut21258804002.png")
> system("convert tmp/8zuez1258804002.ps tmp/8zuez1258804002.png")
> system("convert tmp/99j691258804002.ps tmp/99j691258804002.png")
> system("convert tmp/10nzyh1258804002.ps tmp/10nzyh1258804002.png")
>
>
> proc.time()
user system elapsed
2.509 1.578 2.967