R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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(5.50,235.1,5.40,280.7,5.90,264.6,5.80,240.7,5.10,201.4,4.10,240.8,4.40,241.1,3.60,223.8,3.50,206.1,3.10,174.7,2.90,203.3,2.20,220.5,1.40,299.5,1.20,347.4,1.30,338.3,1.30,327.7,1.30,351.6,1.80,396.6,1.80,438.8,1.80,395.6,1.70,363.5,2.10,378.8,2.00,357.0,1.70,369.0,1.90,464.8,2.30,479.1,2.40,431.3,2.50,366.5,2.80,326.3,2.60,355.1,2.20,331.6,2.80,261.3,2.80,249.0,2.80,205.5,2.30,235.6,2.20,240.9,3.00,264.9,2.90,253.8,2.70,232.3,2.70,193.8,2.30,177.0,2.40,213.2,2.80,207.2,2.30,180.6,2.00,188.6,1.90,175.4,2.30,199.0,2.70,179.6,1.80,225.8,2.00,234.0,2.10,200.2,2.00,183.6,2.40,178.2,1.70,203.2,1.00,208.5,1.20,191.8,1.40,172.8,1.70,148.0,1.80,159.4,1.40,154.5,1.70,213.2,1.60,196.4,1.40,182.8,1.50,176.4,0.90,153.6,1.50,173.2,1.70,171.0,1.60,151.2,1.20,161.9),dim=c(2,69),dimnames=list(c('HIPC','werkloosheid'),1:69))
> y <- array(NA,dim=c(2,69),dimnames=list(c('HIPC','werkloosheid'),1:69))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '2'
> #'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
> 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
werkloosheid HIPC
1 235.1 5.5
2 280.7 5.4
3 264.6 5.9
4 240.7 5.8
5 201.4 5.1
6 240.8 4.1
7 241.1 4.4
8 223.8 3.6
9 206.1 3.5
10 174.7 3.1
11 203.3 2.9
12 220.5 2.2
13 299.5 1.4
14 347.4 1.2
15 338.3 1.3
16 327.7 1.3
17 351.6 1.3
18 396.6 1.8
19 438.8 1.8
20 395.6 1.8
21 363.5 1.7
22 378.8 2.1
23 357.0 2.0
24 369.0 1.7
25 464.8 1.9
26 479.1 2.3
27 431.3 2.4
28 366.5 2.5
29 326.3 2.8
30 355.1 2.6
31 331.6 2.2
32 261.3 2.8
33 249.0 2.8
34 205.5 2.8
35 235.6 2.3
36 240.9 2.2
37 264.9 3.0
38 253.8 2.9
39 232.3 2.7
40 193.8 2.7
41 177.0 2.3
42 213.2 2.4
43 207.2 2.8
44 180.6 2.3
45 188.6 2.0
46 175.4 1.9
47 199.0 2.3
48 179.6 2.7
49 225.8 1.8
50 234.0 2.0
51 200.2 2.1
52 183.6 2.0
53 178.2 2.4
54 203.2 1.7
55 208.5 1.0
56 191.8 1.2
57 172.8 1.4
58 148.0 1.7
59 159.4 1.8
60 154.5 1.4
61 213.2 1.7
62 196.4 1.6
63 182.8 1.4
64 176.4 1.5
65 153.6 0.9
66 173.2 1.5
67 171.0 1.7
68 151.2 1.6
69 161.9 1.2
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) HIPC
254.858 -1.441
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-104.41 -68.38 -26.46 74.72 227.56
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 254.858 24.274 10.499 8.59e-16 ***
HIPC -1.441 9.230 -0.156 0.876
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 86.08 on 67 degrees of freedom
Multiple R-squared: 0.0003637, Adjusted R-squared: -0.01456
F-statistic: 0.02438 on 1 and 67 DF, p-value: 0.8764
> 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.680543e-02 9.361086e-02 9.531946e-01
[2,] 1.969270e-02 3.938540e-02 9.803073e-01
[3,] 5.454464e-03 1.090893e-02 9.945455e-01
[4,] 1.373262e-03 2.746524e-03 9.986267e-01
[5,] 3.993849e-04 7.987697e-04 9.996006e-01
[6,] 2.219532e-04 4.439063e-04 9.997780e-01
[7,] 5.652981e-05 1.130596e-04 9.999435e-01
[8,] 2.915284e-05 5.830568e-05 9.999708e-01
[9,] 7.042853e-04 1.408571e-03 9.992957e-01
[10,] 4.463093e-03 8.926186e-03 9.955369e-01
[11,] 5.384213e-03 1.076843e-02 9.946158e-01
[12,] 4.051588e-03 8.103177e-03 9.959484e-01
[13,] 4.291706e-03 8.583413e-03 9.957083e-01
[14,] 1.274602e-02 2.549204e-02 9.872540e-01
[15,] 6.314638e-02 1.262928e-01 9.368536e-01
[16,] 9.040473e-02 1.808095e-01 9.095953e-01
[17,] 9.030324e-02 1.806065e-01 9.096968e-01
[18,] 1.056771e-01 2.113543e-01 8.943229e-01
[19,] 1.030891e-01 2.061782e-01 8.969109e-01
[20,] 1.211377e-01 2.422753e-01 8.788623e-01
[21,] 4.767259e-01 9.534519e-01 5.232741e-01
[22,] 9.305197e-01 1.389606e-01 6.948031e-02
[23,] 9.960159e-01 7.968129e-03 3.984064e-03
[24,] 9.993374e-01 1.325199e-03 6.625997e-04
[25,] 9.996329e-01 7.341363e-04 3.670682e-04
[26,] 9.999852e-01 2.963207e-05 1.481603e-05
[27,] 9.999999e-01 1.057269e-07 5.286344e-08
[28,] 1.000000e+00 8.879351e-08 4.439675e-08
[29,] 1.000000e+00 9.701873e-08 4.850936e-08
[30,] 1.000000e+00 8.782554e-08 4.391277e-08
[31,] 1.000000e+00 6.774341e-08 3.387170e-08
[32,] 1.000000e+00 3.600821e-08 1.800410e-08
[33,] 1.000000e+00 2.023096e-08 1.011548e-08
[34,] 1.000000e+00 1.010568e-08 5.052838e-09
[35,] 1.000000e+00 8.538773e-09 4.269387e-09
[36,] 1.000000e+00 1.296766e-08 6.483829e-09
[37,] 1.000000e+00 1.201096e-08 6.005478e-09
[38,] 1.000000e+00 1.857819e-08 9.289093e-09
[39,] 1.000000e+00 4.107498e-08 2.053749e-08
[40,] 1.000000e+00 6.488635e-08 3.244317e-08
[41,] 9.999999e-01 1.220824e-07 6.104118e-08
[42,] 9.999999e-01 2.002161e-07 1.001080e-07
[43,] 9.999998e-01 4.833281e-07 2.416640e-07
[44,] 9.999995e-01 9.063337e-07 4.531669e-07
[45,] 9.999996e-01 7.354241e-07 3.677120e-07
[46,] 9.999999e-01 2.381061e-07 1.190531e-07
[47,] 9.999998e-01 4.971321e-07 2.485660e-07
[48,] 9.999993e-01 1.460969e-06 7.304847e-07
[49,] 9.999977e-01 4.615425e-06 2.307713e-06
[50,] 9.999966e-01 6.833652e-06 3.416826e-06
[51,] 9.999971e-01 5.882217e-06 2.941109e-06
[52,] 9.999949e-01 1.016029e-05 5.080144e-06
[53,] 9.999815e-01 3.708029e-05 1.854014e-05
[54,] 9.999748e-01 5.046763e-05 2.523382e-05
[55,] 9.999569e-01 8.627448e-05 4.313724e-05
[56,] 9.998784e-01 2.431490e-04 1.215745e-04
[57,] 9.998846e-01 2.307426e-04 1.153713e-04
[58,] 9.998254e-01 3.491474e-04 1.745737e-04
[59,] 9.995305e-01 9.390129e-04 4.695065e-04
[60,] 9.976204e-01 4.759122e-03 2.379561e-03
> postscript(file="/var/wessaorg/rcomp/tmp/19a4b1324307058.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/2gln21324307058.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/3lco61324307058.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/4kw0u1324307058.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/5gape1324307058.ps",horizontal=F,onefile=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
-11.832011 33.623889 18.244388 -5.799712 -46.108411 -8.149410
7 8 9 10 11 12
-7.417110 -25.869910 -43.714010 -75.690409 -47.378609 -31.187309
13 14 15 16 17 18
46.659892 94.271692 85.315792 74.715792 98.615792 144.336292
19 20 21 22 23 24
186.536292 143.336292 111.092192 126.968592 105.024492 116.592192
25 26 27 28 29 30
212.680392 227.556791 179.900891 115.244991 75.477291 103.989091
31 32 33 34 35 36
79.912691 10.477291 -1.822709 -45.322709 -15.943209 -10.787309
37 38 39 40 41 42
14.365491 3.121391 -18.666809 -57.166809 -74.543209 -38.199109
43 44 45 46 47 48
-43.622709 -70.943209 -63.375508 -76.719608 -52.543209 -71.366809
49 50 51 52 53 54
-26.463708 -17.975508 -51.631408 -68.375508 -73.199109 -49.207808
55 56 57 58 59 60
-44.916508 -61.328308 -80.040108 -104.407808 -92.863708 -98.340108
61 62 63 64 65 66
-39.207808 -56.151908 -70.040108 -76.296008 -99.960607 -79.496008
67 68 69
-81.407808 -101.351908 -91.228308
> postscript(file="/var/wessaorg/rcomp/tmp/68mry1324307058.ps",horizontal=F,onefile=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 -11.832011 NA
1 33.623889 -11.832011
2 18.244388 33.623889
3 -5.799712 18.244388
4 -46.108411 -5.799712
5 -8.149410 -46.108411
6 -7.417110 -8.149410
7 -25.869910 -7.417110
8 -43.714010 -25.869910
9 -75.690409 -43.714010
10 -47.378609 -75.690409
11 -31.187309 -47.378609
12 46.659892 -31.187309
13 94.271692 46.659892
14 85.315792 94.271692
15 74.715792 85.315792
16 98.615792 74.715792
17 144.336292 98.615792
18 186.536292 144.336292
19 143.336292 186.536292
20 111.092192 143.336292
21 126.968592 111.092192
22 105.024492 126.968592
23 116.592192 105.024492
24 212.680392 116.592192
25 227.556791 212.680392
26 179.900891 227.556791
27 115.244991 179.900891
28 75.477291 115.244991
29 103.989091 75.477291
30 79.912691 103.989091
31 10.477291 79.912691
32 -1.822709 10.477291
33 -45.322709 -1.822709
34 -15.943209 -45.322709
35 -10.787309 -15.943209
36 14.365491 -10.787309
37 3.121391 14.365491
38 -18.666809 3.121391
39 -57.166809 -18.666809
40 -74.543209 -57.166809
41 -38.199109 -74.543209
42 -43.622709 -38.199109
43 -70.943209 -43.622709
44 -63.375508 -70.943209
45 -76.719608 -63.375508
46 -52.543209 -76.719608
47 -71.366809 -52.543209
48 -26.463708 -71.366809
49 -17.975508 -26.463708
50 -51.631408 -17.975508
51 -68.375508 -51.631408
52 -73.199109 -68.375508
53 -49.207808 -73.199109
54 -44.916508 -49.207808
55 -61.328308 -44.916508
56 -80.040108 -61.328308
57 -104.407808 -80.040108
58 -92.863708 -104.407808
59 -98.340108 -92.863708
60 -39.207808 -98.340108
61 -56.151908 -39.207808
62 -70.040108 -56.151908
63 -76.296008 -70.040108
64 -99.960607 -76.296008
65 -79.496008 -99.960607
66 -81.407808 -79.496008
67 -101.351908 -81.407808
68 -91.228308 -101.351908
69 NA -91.228308
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 33.623889 -11.832011
[2,] 18.244388 33.623889
[3,] -5.799712 18.244388
[4,] -46.108411 -5.799712
[5,] -8.149410 -46.108411
[6,] -7.417110 -8.149410
[7,] -25.869910 -7.417110
[8,] -43.714010 -25.869910
[9,] -75.690409 -43.714010
[10,] -47.378609 -75.690409
[11,] -31.187309 -47.378609
[12,] 46.659892 -31.187309
[13,] 94.271692 46.659892
[14,] 85.315792 94.271692
[15,] 74.715792 85.315792
[16,] 98.615792 74.715792
[17,] 144.336292 98.615792
[18,] 186.536292 144.336292
[19,] 143.336292 186.536292
[20,] 111.092192 143.336292
[21,] 126.968592 111.092192
[22,] 105.024492 126.968592
[23,] 116.592192 105.024492
[24,] 212.680392 116.592192
[25,] 227.556791 212.680392
[26,] 179.900891 227.556791
[27,] 115.244991 179.900891
[28,] 75.477291 115.244991
[29,] 103.989091 75.477291
[30,] 79.912691 103.989091
[31,] 10.477291 79.912691
[32,] -1.822709 10.477291
[33,] -45.322709 -1.822709
[34,] -15.943209 -45.322709
[35,] -10.787309 -15.943209
[36,] 14.365491 -10.787309
[37,] 3.121391 14.365491
[38,] -18.666809 3.121391
[39,] -57.166809 -18.666809
[40,] -74.543209 -57.166809
[41,] -38.199109 -74.543209
[42,] -43.622709 -38.199109
[43,] -70.943209 -43.622709
[44,] -63.375508 -70.943209
[45,] -76.719608 -63.375508
[46,] -52.543209 -76.719608
[47,] -71.366809 -52.543209
[48,] -26.463708 -71.366809
[49,] -17.975508 -26.463708
[50,] -51.631408 -17.975508
[51,] -68.375508 -51.631408
[52,] -73.199109 -68.375508
[53,] -49.207808 -73.199109
[54,] -44.916508 -49.207808
[55,] -61.328308 -44.916508
[56,] -80.040108 -61.328308
[57,] -104.407808 -80.040108
[58,] -92.863708 -104.407808
[59,] -98.340108 -92.863708
[60,] -39.207808 -98.340108
[61,] -56.151908 -39.207808
[62,] -70.040108 -56.151908
[63,] -76.296008 -70.040108
[64,] -99.960607 -76.296008
[65,] -79.496008 -99.960607
[66,] -81.407808 -79.496008
[67,] -101.351908 -81.407808
[68,] -91.228308 -101.351908
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 33.623889 -11.832011
2 18.244388 33.623889
3 -5.799712 18.244388
4 -46.108411 -5.799712
5 -8.149410 -46.108411
6 -7.417110 -8.149410
7 -25.869910 -7.417110
8 -43.714010 -25.869910
9 -75.690409 -43.714010
10 -47.378609 -75.690409
11 -31.187309 -47.378609
12 46.659892 -31.187309
13 94.271692 46.659892
14 85.315792 94.271692
15 74.715792 85.315792
16 98.615792 74.715792
17 144.336292 98.615792
18 186.536292 144.336292
19 143.336292 186.536292
20 111.092192 143.336292
21 126.968592 111.092192
22 105.024492 126.968592
23 116.592192 105.024492
24 212.680392 116.592192
25 227.556791 212.680392
26 179.900891 227.556791
27 115.244991 179.900891
28 75.477291 115.244991
29 103.989091 75.477291
30 79.912691 103.989091
31 10.477291 79.912691
32 -1.822709 10.477291
33 -45.322709 -1.822709
34 -15.943209 -45.322709
35 -10.787309 -15.943209
36 14.365491 -10.787309
37 3.121391 14.365491
38 -18.666809 3.121391
39 -57.166809 -18.666809
40 -74.543209 -57.166809
41 -38.199109 -74.543209
42 -43.622709 -38.199109
43 -70.943209 -43.622709
44 -63.375508 -70.943209
45 -76.719608 -63.375508
46 -52.543209 -76.719608
47 -71.366809 -52.543209
48 -26.463708 -71.366809
49 -17.975508 -26.463708
50 -51.631408 -17.975508
51 -68.375508 -51.631408
52 -73.199109 -68.375508
53 -49.207808 -73.199109
54 -44.916508 -49.207808
55 -61.328308 -44.916508
56 -80.040108 -61.328308
57 -104.407808 -80.040108
58 -92.863708 -104.407808
59 -98.340108 -92.863708
60 -39.207808 -98.340108
61 -56.151908 -39.207808
62 -70.040108 -56.151908
63 -76.296008 -70.040108
64 -99.960607 -76.296008
65 -79.496008 -99.960607
66 -81.407808 -79.496008
67 -101.351908 -81.407808
68 -91.228308 -101.351908
> 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/wessaorg/rcomp/tmp/7l7do1324307058.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/8q6l01324307058.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/9gpbx1324307058.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/10ayzw1324307058.ps",horizontal=F,onefile=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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11bo2v1324307058.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/wessaorg/rcomp/tmp/12k3q11324307058.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/wessaorg/rcomp/tmp/13dffl1324307058.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/wessaorg/rcomp/tmp/1426kn1324307058.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/wessaorg/rcomp/tmp/15dyc41324307058.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/wessaorg/rcomp/tmp/164y4k1324307058.tab")
+ }
>
> try(system("convert tmp/19a4b1324307058.ps tmp/19a4b1324307058.png",intern=TRUE))
character(0)
> try(system("convert tmp/2gln21324307058.ps tmp/2gln21324307058.png",intern=TRUE))
character(0)
> try(system("convert tmp/3lco61324307058.ps tmp/3lco61324307058.png",intern=TRUE))
character(0)
> try(system("convert tmp/4kw0u1324307058.ps tmp/4kw0u1324307058.png",intern=TRUE))
character(0)
> try(system("convert tmp/5gape1324307058.ps tmp/5gape1324307058.png",intern=TRUE))
character(0)
> try(system("convert tmp/68mry1324307058.ps tmp/68mry1324307058.png",intern=TRUE))
character(0)
> try(system("convert tmp/7l7do1324307058.ps tmp/7l7do1324307058.png",intern=TRUE))
character(0)
> try(system("convert tmp/8q6l01324307058.ps tmp/8q6l01324307058.png",intern=TRUE))
character(0)
> try(system("convert tmp/9gpbx1324307058.ps tmp/9gpbx1324307058.png",intern=TRUE))
character(0)
> try(system("convert tmp/10ayzw1324307058.ps tmp/10ayzw1324307058.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.278 0.663 3.956