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(593530.00,0,610943.00,0,612613.00,0,611324.00,0,594167.00,0,595454.00,0,590865.00,0,589379.00,0,584428.00,0,573100.00,0,567456.00,0,569028.00,0,620735.00,0,628884.00,0,628232.00,0,612117.00,0,595404.00,0,597141.00,0,593408.00,0,590072.00,0,579799.00,0,574205.00,0,572775.00,0,572942.00,0,619567.00,0,625809.00,0,619916.00,0,587625.00,0,565742.00,0,557274.00,0,560576.00,0,548854.00,0,531673.00,0,525919.00,0,511038.00,0,498662.00,0,555362.00,0,564591.00,0,541657.00,0,527070.00,0,509846.00,0,514258.00,0,516922.00,0,507561.00,0,492622.00,0,490243.00,0,469357.00,0,477580.00,0,528379.00,1,533590.00,1,517945.00,1,506174.00,1,501866.00,1,516141.00,1,528222.00,1,532638.00,1,536322.00,1,536535.00,1,523597.00,1,536214.00,1,586570.00,1,596594.00,1,580523.00,1),dim=c(2,63),dimnames=list(c('werklozen','crisis
'),1:63))
> y <- array(NA,dim=c(2,63),dimnames=list(c('werklozen','crisis
'),1:63))
> 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
werklozen crisis\r M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 593530 0 1 0 0 0 0 0 0 0 0 0 0
2 610943 0 0 1 0 0 0 0 0 0 0 0 0
3 612613 0 0 0 1 0 0 0 0 0 0 0 0
4 611324 0 0 0 0 1 0 0 0 0 0 0 0
5 594167 0 0 0 0 0 1 0 0 0 0 0 0
6 595454 0 0 0 0 0 0 1 0 0 0 0 0
7 590865 0 0 0 0 0 0 0 1 0 0 0 0
8 589379 0 0 0 0 0 0 0 0 1 0 0 0
9 584428 0 0 0 0 0 0 0 0 0 1 0 0
10 573100 0 0 0 0 0 0 0 0 0 0 1 0
11 567456 0 0 0 0 0 0 0 0 0 0 0 1
12 569028 0 0 0 0 0 0 0 0 0 0 0 0
13 620735 0 1 0 0 0 0 0 0 0 0 0 0
14 628884 0 0 1 0 0 0 0 0 0 0 0 0
15 628232 0 0 0 1 0 0 0 0 0 0 0 0
16 612117 0 0 0 0 1 0 0 0 0 0 0 0
17 595404 0 0 0 0 0 1 0 0 0 0 0 0
18 597141 0 0 0 0 0 0 1 0 0 0 0 0
19 593408 0 0 0 0 0 0 0 1 0 0 0 0
20 590072 0 0 0 0 0 0 0 0 1 0 0 0
21 579799 0 0 0 0 0 0 0 0 0 1 0 0
22 574205 0 0 0 0 0 0 0 0 0 0 1 0
23 572775 0 0 0 0 0 0 0 0 0 0 0 1
24 572942 0 0 0 0 0 0 0 0 0 0 0 0
25 619567 0 1 0 0 0 0 0 0 0 0 0 0
26 625809 0 0 1 0 0 0 0 0 0 0 0 0
27 619916 0 0 0 1 0 0 0 0 0 0 0 0
28 587625 0 0 0 0 1 0 0 0 0 0 0 0
29 565742 0 0 0 0 0 1 0 0 0 0 0 0
30 557274 0 0 0 0 0 0 1 0 0 0 0 0
31 560576 0 0 0 0 0 0 0 1 0 0 0 0
32 548854 0 0 0 0 0 0 0 0 1 0 0 0
33 531673 0 0 0 0 0 0 0 0 0 1 0 0
34 525919 0 0 0 0 0 0 0 0 0 0 1 0
35 511038 0 0 0 0 0 0 0 0 0 0 0 1
36 498662 0 0 0 0 0 0 0 0 0 0 0 0
37 555362 0 1 0 0 0 0 0 0 0 0 0 0
38 564591 0 0 1 0 0 0 0 0 0 0 0 0
39 541657 0 0 0 1 0 0 0 0 0 0 0 0
40 527070 0 0 0 0 1 0 0 0 0 0 0 0
41 509846 0 0 0 0 0 1 0 0 0 0 0 0
42 514258 0 0 0 0 0 0 1 0 0 0 0 0
43 516922 0 0 0 0 0 0 0 1 0 0 0 0
44 507561 0 0 0 0 0 0 0 0 1 0 0 0
45 492622 0 0 0 0 0 0 0 0 0 1 0 0
46 490243 0 0 0 0 0 0 0 0 0 0 1 0
47 469357 0 0 0 0 0 0 0 0 0 0 0 1
48 477580 0 0 0 0 0 0 0 0 0 0 0 0
49 528379 1 1 0 0 0 0 0 0 0 0 0 0
50 533590 1 0 1 0 0 0 0 0 0 0 0 0
51 517945 1 0 0 1 0 0 0 0 0 0 0 0
52 506174 1 0 0 0 1 0 0 0 0 0 0 0
53 501866 1 0 0 0 0 1 0 0 0 0 0 0
54 516141 1 0 0 0 0 0 1 0 0 0 0 0
55 528222 1 0 0 0 0 0 0 1 0 0 0 0
56 532638 1 0 0 0 0 0 0 0 1 0 0 0
57 536322 1 0 0 0 0 0 0 0 0 1 0 0
58 536535 1 0 0 0 0 0 0 0 0 0 1 0
59 523597 1 0 0 0 0 0 0 0 0 0 0 1
60 536214 1 0 0 0 0 0 0 0 0 0 0 0
61 586570 1 1 0 0 0 0 0 0 0 0 0 0
62 596594 1 0 1 0 0 0 0 0 0 0 0 0
63 580523 1 0 0 1 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) `crisis\r` M1 M2 M3 M4
537943 -35288 57844 67222 57301 37977
M5 M6 M7 M8 M9 M10
22520 25168 27113 22816 14084 9115
M11
-2041
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-66545 -33289 17369 28062 36873
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 537943 16836 31.952 < 2e-16 ***
`crisis\r` -35288 11150 -3.165 0.00264 **
M1 57844 22644 2.555 0.01373 *
M2 67222 22644 2.969 0.00458 **
M3 57301 22644 2.531 0.01458 *
M4 37977 23600 1.609 0.11387
M5 22520 23600 0.954 0.34455
M6 25168 23600 1.066 0.29133
M7 27113 23600 1.149 0.25607
M8 22816 23600 0.967 0.33831
M9 14084 23600 0.597 0.55335
M10 9115 23600 0.386 0.70095
M11 -2041 23600 -0.086 0.93144
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 37310 on 50 degrees of freedom
Multiple R-squared: 0.3679, Adjusted R-squared: 0.2161
F-statistic: 2.425 on 12 and 50 DF, p-value: 0.01440
> 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,] 6.856811e-02 1.371362e-01 0.9314319
[2,] 2.298249e-02 4.596498e-02 0.9770175
[3,] 7.360790e-03 1.472158e-02 0.9926392
[4,] 2.203895e-03 4.407791e-03 0.9977961
[5,] 6.435883e-04 1.287177e-03 0.9993564
[6,] 1.964215e-04 3.928431e-04 0.9998036
[7,] 5.567986e-05 1.113597e-04 0.9999443
[8,] 2.109207e-05 4.218415e-05 0.9999789
[9,] 8.167038e-06 1.633408e-05 0.9999918
[10,] 6.254677e-06 1.250935e-05 0.9999937
[11,] 3.049239e-06 6.098477e-06 0.9999970
[12,] 2.214457e-06 4.428914e-06 0.9999978
[13,] 2.464122e-05 4.928243e-05 0.9999754
[14,] 2.831895e-04 5.663789e-04 0.9997168
[15,] 3.351688e-03 6.703377e-03 0.9966483
[16,] 9.369059e-03 1.873812e-02 0.9906309
[17,] 2.984397e-02 5.968794e-02 0.9701560
[18,] 8.257297e-02 1.651459e-01 0.9174270
[19,] 1.334877e-01 2.669753e-01 0.8665123
[20,] 2.360969e-01 4.721938e-01 0.7639031
[21,] 3.673243e-01 7.346486e-01 0.6326757
[22,] 4.130432e-01 8.260864e-01 0.5869568
[23,] 4.503870e-01 9.007741e-01 0.5496130
[24,] 5.350296e-01 9.299408e-01 0.4649704
[25,] 6.520387e-01 6.959227e-01 0.3479613
[26,] 7.074577e-01 5.850846e-01 0.2925423
[27,] 7.180023e-01 5.639954e-01 0.2819977
[28,] 6.983281e-01 6.033437e-01 0.3016719
[29,] 6.489506e-01 7.020987e-01 0.3510494
[30,] 5.635838e-01 8.728324e-01 0.4364162
[31,] 4.508941e-01 9.017881e-01 0.5491059
[32,] 3.325623e-01 6.651246e-01 0.6674377
> postscript(file="/var/www/html/rcomp/tmp/1te141258646290.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/2q5vj1258646290.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/3nvpm1258646290.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/40k1c1258646290.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/5ifsm1258646290.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 = 63
Frequency = 1
1 2 3 4 5 6 7
-2256.403 5778.597 17369.431 35404.458 33704.458 32342.858 25808.858
8 9 10 11 12 13 14
28620.658 32401.658 26042.058 31553.858 31085.258 24948.597 23719.597
15 16 17 18 19 20 21
32988.431 36197.458 34941.458 34029.858 28351.858 29313.658 27772.658
22 23 24 25 26 27 28
27147.058 36872.858 34999.258 23780.597 20644.597 24672.431 11705.458
29 30 31 32 33 34 35
5279.458 -5837.142 -4480.142 -11904.342 -20353.342 -21138.942 -24864.142
36 37 38 39 40 41 42
-39280.742 -40424.403 -40573.403 -53586.569 -48849.542 -50616.542 -48853.142
43 44 45 46 47 48 49
-48134.142 -53197.342 -59404.342 -56814.942 -66545.142 -60362.742 -32119.694
50 51 52 53 54 55 56
-36286.694 -42010.861 -34457.833 -23308.833 -11682.433 -1546.433 7167.367
57 58 59 60 61 62 63
19583.367 24764.767 22982.567 33558.967 26071.306 26717.306 20567.139
> postscript(file="/var/www/html/rcomp/tmp/6ydbo1258646290.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 = 63
Frequency = 1
lag(myerror, k = 1) myerror
0 -2256.403 NA
1 5778.597 -2256.403
2 17369.431 5778.597
3 35404.458 17369.431
4 33704.458 35404.458
5 32342.858 33704.458
6 25808.858 32342.858
7 28620.658 25808.858
8 32401.658 28620.658
9 26042.058 32401.658
10 31553.858 26042.058
11 31085.258 31553.858
12 24948.597 31085.258
13 23719.597 24948.597
14 32988.431 23719.597
15 36197.458 32988.431
16 34941.458 36197.458
17 34029.858 34941.458
18 28351.858 34029.858
19 29313.658 28351.858
20 27772.658 29313.658
21 27147.058 27772.658
22 36872.858 27147.058
23 34999.258 36872.858
24 23780.597 34999.258
25 20644.597 23780.597
26 24672.431 20644.597
27 11705.458 24672.431
28 5279.458 11705.458
29 -5837.142 5279.458
30 -4480.142 -5837.142
31 -11904.342 -4480.142
32 -20353.342 -11904.342
33 -21138.942 -20353.342
34 -24864.142 -21138.942
35 -39280.742 -24864.142
36 -40424.403 -39280.742
37 -40573.403 -40424.403
38 -53586.569 -40573.403
39 -48849.542 -53586.569
40 -50616.542 -48849.542
41 -48853.142 -50616.542
42 -48134.142 -48853.142
43 -53197.342 -48134.142
44 -59404.342 -53197.342
45 -56814.942 -59404.342
46 -66545.142 -56814.942
47 -60362.742 -66545.142
48 -32119.694 -60362.742
49 -36286.694 -32119.694
50 -42010.861 -36286.694
51 -34457.833 -42010.861
52 -23308.833 -34457.833
53 -11682.433 -23308.833
54 -1546.433 -11682.433
55 7167.367 -1546.433
56 19583.367 7167.367
57 24764.767 19583.367
58 22982.567 24764.767
59 33558.967 22982.567
60 26071.306 33558.967
61 26717.306 26071.306
62 20567.139 26717.306
63 NA 20567.139
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 5778.597 -2256.403
[2,] 17369.431 5778.597
[3,] 35404.458 17369.431
[4,] 33704.458 35404.458
[5,] 32342.858 33704.458
[6,] 25808.858 32342.858
[7,] 28620.658 25808.858
[8,] 32401.658 28620.658
[9,] 26042.058 32401.658
[10,] 31553.858 26042.058
[11,] 31085.258 31553.858
[12,] 24948.597 31085.258
[13,] 23719.597 24948.597
[14,] 32988.431 23719.597
[15,] 36197.458 32988.431
[16,] 34941.458 36197.458
[17,] 34029.858 34941.458
[18,] 28351.858 34029.858
[19,] 29313.658 28351.858
[20,] 27772.658 29313.658
[21,] 27147.058 27772.658
[22,] 36872.858 27147.058
[23,] 34999.258 36872.858
[24,] 23780.597 34999.258
[25,] 20644.597 23780.597
[26,] 24672.431 20644.597
[27,] 11705.458 24672.431
[28,] 5279.458 11705.458
[29,] -5837.142 5279.458
[30,] -4480.142 -5837.142
[31,] -11904.342 -4480.142
[32,] -20353.342 -11904.342
[33,] -21138.942 -20353.342
[34,] -24864.142 -21138.942
[35,] -39280.742 -24864.142
[36,] -40424.403 -39280.742
[37,] -40573.403 -40424.403
[38,] -53586.569 -40573.403
[39,] -48849.542 -53586.569
[40,] -50616.542 -48849.542
[41,] -48853.142 -50616.542
[42,] -48134.142 -48853.142
[43,] -53197.342 -48134.142
[44,] -59404.342 -53197.342
[45,] -56814.942 -59404.342
[46,] -66545.142 -56814.942
[47,] -60362.742 -66545.142
[48,] -32119.694 -60362.742
[49,] -36286.694 -32119.694
[50,] -42010.861 -36286.694
[51,] -34457.833 -42010.861
[52,] -23308.833 -34457.833
[53,] -11682.433 -23308.833
[54,] -1546.433 -11682.433
[55,] 7167.367 -1546.433
[56,] 19583.367 7167.367
[57,] 24764.767 19583.367
[58,] 22982.567 24764.767
[59,] 33558.967 22982.567
[60,] 26071.306 33558.967
[61,] 26717.306 26071.306
[62,] 20567.139 26717.306
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 5778.597 -2256.403
2 17369.431 5778.597
3 35404.458 17369.431
4 33704.458 35404.458
5 32342.858 33704.458
6 25808.858 32342.858
7 28620.658 25808.858
8 32401.658 28620.658
9 26042.058 32401.658
10 31553.858 26042.058
11 31085.258 31553.858
12 24948.597 31085.258
13 23719.597 24948.597
14 32988.431 23719.597
15 36197.458 32988.431
16 34941.458 36197.458
17 34029.858 34941.458
18 28351.858 34029.858
19 29313.658 28351.858
20 27772.658 29313.658
21 27147.058 27772.658
22 36872.858 27147.058
23 34999.258 36872.858
24 23780.597 34999.258
25 20644.597 23780.597
26 24672.431 20644.597
27 11705.458 24672.431
28 5279.458 11705.458
29 -5837.142 5279.458
30 -4480.142 -5837.142
31 -11904.342 -4480.142
32 -20353.342 -11904.342
33 -21138.942 -20353.342
34 -24864.142 -21138.942
35 -39280.742 -24864.142
36 -40424.403 -39280.742
37 -40573.403 -40424.403
38 -53586.569 -40573.403
39 -48849.542 -53586.569
40 -50616.542 -48849.542
41 -48853.142 -50616.542
42 -48134.142 -48853.142
43 -53197.342 -48134.142
44 -59404.342 -53197.342
45 -56814.942 -59404.342
46 -66545.142 -56814.942
47 -60362.742 -66545.142
48 -32119.694 -60362.742
49 -36286.694 -32119.694
50 -42010.861 -36286.694
51 -34457.833 -42010.861
52 -23308.833 -34457.833
53 -11682.433 -23308.833
54 -1546.433 -11682.433
55 7167.367 -1546.433
56 19583.367 7167.367
57 24764.767 19583.367
58 22982.567 24764.767
59 33558.967 22982.567
60 26071.306 33558.967
61 26717.306 26071.306
62 20567.139 26717.306
> 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/7ak711258646290.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/8ah571258646290.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/90vff1258646290.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/10xm0z1258646290.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/11lb5q1258646290.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/12ha8s1258646290.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/13nm5a1258646290.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/14kwv21258646290.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/15xzko1258646290.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/16q2p11258646290.tab")
+ }
>
> system("convert tmp/1te141258646290.ps tmp/1te141258646290.png")
> system("convert tmp/2q5vj1258646290.ps tmp/2q5vj1258646290.png")
> system("convert tmp/3nvpm1258646290.ps tmp/3nvpm1258646290.png")
> system("convert tmp/40k1c1258646290.ps tmp/40k1c1258646290.png")
> system("convert tmp/5ifsm1258646290.ps tmp/5ifsm1258646290.png")
> system("convert tmp/6ydbo1258646290.ps tmp/6ydbo1258646290.png")
> system("convert tmp/7ak711258646290.ps tmp/7ak711258646290.png")
> system("convert tmp/8ah571258646290.ps tmp/8ah571258646290.png")
> system("convert tmp/90vff1258646290.ps tmp/90vff1258646290.png")
> system("convert tmp/10xm0z1258646290.ps tmp/10xm0z1258646290.png")
>
>
> proc.time()
user system elapsed
2.437 1.563 2.849