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(441,1919,449,1911,452,1870,462,2263,455,1802,461,1863,461,1989,463,2197,462,2409,456,2502,455,2593,456,2598,472,2053,472,2213,471,2238,465,2359,459,2151,465,2474,468,3079,467,2312,463,2565,460,1972,462,2484,461,2202,476,2151,476,1976,471,2012,453,2114,443,1772,442,1957,444,2070,438,1990,427,2182,424,2008,416,1916,406,2397,431,2114,434,1778,418,1641,412,2186,404,1773,409,1785,412,2217,406,2153,398,1895,397,2475,385,1793,390,2308,413,2051,413,1898,401,2142,397,1874,397,1560,409,1808,419,1575,424,1525,428,1997,430,1753,424,1623,433,2251,456,1890),dim=c(2,61),dimnames=list(c('wkl','bvg'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('wkl','bvg'),1:61))
> 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
wkl bvg M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 441 1919 1 0 0 0 0 0 0 0 0 0 0 1
2 449 1911 0 1 0 0 0 0 0 0 0 0 0 2
3 452 1870 0 0 1 0 0 0 0 0 0 0 0 3
4 462 2263 0 0 0 1 0 0 0 0 0 0 0 4
5 455 1802 0 0 0 0 1 0 0 0 0 0 0 5
6 461 1863 0 0 0 0 0 1 0 0 0 0 0 6
7 461 1989 0 0 0 0 0 0 1 0 0 0 0 7
8 463 2197 0 0 0 0 0 0 0 1 0 0 0 8
9 462 2409 0 0 0 0 0 0 0 0 1 0 0 9
10 456 2502 0 0 0 0 0 0 0 0 0 1 0 10
11 455 2593 0 0 0 0 0 0 0 0 0 0 1 11
12 456 2598 0 0 0 0 0 0 0 0 0 0 0 12
13 472 2053 1 0 0 0 0 0 0 0 0 0 0 13
14 472 2213 0 1 0 0 0 0 0 0 0 0 0 14
15 471 2238 0 0 1 0 0 0 0 0 0 0 0 15
16 465 2359 0 0 0 1 0 0 0 0 0 0 0 16
17 459 2151 0 0 0 0 1 0 0 0 0 0 0 17
18 465 2474 0 0 0 0 0 1 0 0 0 0 0 18
19 468 3079 0 0 0 0 0 0 1 0 0 0 0 19
20 467 2312 0 0 0 0 0 0 0 1 0 0 0 20
21 463 2565 0 0 0 0 0 0 0 0 1 0 0 21
22 460 1972 0 0 0 0 0 0 0 0 0 1 0 22
23 462 2484 0 0 0 0 0 0 0 0 0 0 1 23
24 461 2202 0 0 0 0 0 0 0 0 0 0 0 24
25 476 2151 1 0 0 0 0 0 0 0 0 0 0 25
26 476 1976 0 1 0 0 0 0 0 0 0 0 0 26
27 471 2012 0 0 1 0 0 0 0 0 0 0 0 27
28 453 2114 0 0 0 1 0 0 0 0 0 0 0 28
29 443 1772 0 0 0 0 1 0 0 0 0 0 0 29
30 442 1957 0 0 0 0 0 1 0 0 0 0 0 30
31 444 2070 0 0 0 0 0 0 1 0 0 0 0 31
32 438 1990 0 0 0 0 0 0 0 1 0 0 0 32
33 427 2182 0 0 0 0 0 0 0 0 1 0 0 33
34 424 2008 0 0 0 0 0 0 0 0 0 1 0 34
35 416 1916 0 0 0 0 0 0 0 0 0 0 1 35
36 406 2397 0 0 0 0 0 0 0 0 0 0 0 36
37 431 2114 1 0 0 0 0 0 0 0 0 0 0 37
38 434 1778 0 1 0 0 0 0 0 0 0 0 0 38
39 418 1641 0 0 1 0 0 0 0 0 0 0 0 39
40 412 2186 0 0 0 1 0 0 0 0 0 0 0 40
41 404 1773 0 0 0 0 1 0 0 0 0 0 0 41
42 409 1785 0 0 0 0 0 1 0 0 0 0 0 42
43 412 2217 0 0 0 0 0 0 1 0 0 0 0 43
44 406 2153 0 0 0 0 0 0 0 1 0 0 0 44
45 398 1895 0 0 0 0 0 0 0 0 1 0 0 45
46 397 2475 0 0 0 0 0 0 0 0 0 1 0 46
47 385 1793 0 0 0 0 0 0 0 0 0 0 1 47
48 390 2308 0 0 0 0 0 0 0 0 0 0 0 48
49 413 2051 1 0 0 0 0 0 0 0 0 0 0 49
50 413 1898 0 1 0 0 0 0 0 0 0 0 0 50
51 401 2142 0 0 1 0 0 0 0 0 0 0 0 51
52 397 1874 0 0 0 1 0 0 0 0 0 0 0 52
53 397 1560 0 0 0 0 1 0 0 0 0 0 0 53
54 409 1808 0 0 0 0 0 1 0 0 0 0 0 54
55 419 1575 0 0 0 0 0 0 1 0 0 0 0 55
56 424 1525 0 0 0 0 0 0 0 1 0 0 0 56
57 428 1997 0 0 0 0 0 0 0 0 1 0 0 57
58 430 1753 0 0 0 0 0 0 0 0 0 1 0 58
59 424 1623 0 0 0 0 0 0 0 0 0 0 1 59
60 433 2251 0 0 0 0 0 0 0 0 0 0 0 60
61 456 1890 1 0 0 0 0 0 0 0 0 0 0 61
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) bvg M1 M2 M3 M4
433.12004 0.01311 18.35638 15.14046 9.57257 3.39628
M5 M6 M7 M8 M9 M10
2.71830 7.10981 8.94023 10.67964 5.36104 5.01235
M11 t
1.76665 -0.96510
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-34.6688 -12.4817 -0.3231 12.1927 38.6172
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 433.12004 30.31547 14.287 < 2e-16 ***
bvg 0.01311 0.01087 1.206 0.234
M1 18.35638 12.18972 1.506 0.139
M2 15.14046 13.16820 1.150 0.256
M3 9.57257 13.01601 0.735 0.466
M4 3.39628 12.39827 0.274 0.785
M5 2.71830 13.69480 0.198 0.844
M6 7.10981 12.90086 0.551 0.584
M7 8.94023 12.25559 0.729 0.469
M8 10.67964 12.62265 0.846 0.402
M9 5.36104 12.17207 0.440 0.662
M10 5.01235 12.28060 0.408 0.685
M11 1.76665 12.40389 0.142 0.887
t -0.96510 0.16470 -5.860 4.38e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 19.02 on 47 degrees of freedom
Multiple R-squared: 0.5858, Adjusted R-squared: 0.4713
F-statistic: 5.114 on 13 and 47 DF, p-value: 1.533e-05
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 1.244905e-01 0.2489810794 0.87550946
[2,] 5.220701e-02 0.1044140237 0.94779299
[3,] 2.067205e-02 0.0413441031 0.97932795
[4,] 1.117375e-02 0.0223475079 0.98882625
[5,] 6.523141e-03 0.0130462811 0.99347686
[6,] 2.675260e-03 0.0053505195 0.99732474
[7,] 1.409291e-03 0.0028185811 0.99859071
[8,] 4.996949e-04 0.0009993898 0.99950031
[9,] 1.991032e-04 0.0003982064 0.99980090
[10,] 9.784247e-05 0.0001956849 0.99990216
[11,] 9.917335e-05 0.0001983467 0.99990083
[12,] 9.538540e-04 0.0019077079 0.99904615
[13,] 4.211127e-03 0.0084222541 0.99578887
[14,] 1.608391e-02 0.0321678268 0.98391609
[15,] 2.575147e-02 0.0515029469 0.97424853
[16,] 5.430613e-02 0.1086122656 0.94569387
[17,] 1.085628e-01 0.2171256210 0.89143719
[18,] 1.315947e-01 0.2631894262 0.86840529
[19,] 1.654019e-01 0.3308037899 0.83459811
[20,] 3.020020e-01 0.6040039312 0.69799803
[21,] 2.798829e-01 0.5597657843 0.72011711
[22,] 2.985596e-01 0.5971192946 0.70144035
[23,] 3.683512e-01 0.7367023733 0.63164881
[24,] 5.591687e-01 0.8816626658 0.44083133
[25,] 7.206222e-01 0.5587555389 0.27937777
[26,] 9.263331e-01 0.1473337285 0.07366686
[27,] 9.526758e-01 0.0946483604 0.04732418
[28,] 9.664082e-01 0.0671835690 0.03359178
> postscript(file="/var/www/html/rcomp/tmp/1p52h1260820143.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/212m61260820143.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/38oln1260820143.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/4qqab1260820143.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/55yda1260820143.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 = 61
Frequency = 1
1 2 3 4 5 6
-34.6687620 -22.3828712 -12.3123897 -0.3231012 0.3635345 2.1374270
7 8 9 10 11 12
-0.3797150 -1.8808396 0.6236152 -5.2817965 -3.2639773 0.4022172
13 14 15 16 17 18
6.1557025 8.2391701 13.4444139 12.9995305 11.3694217 9.7085829
19 20 21 22 23 24
3.9119130 12.1927084 11.1596671 17.2474833 16.7461349 22.1748024
25 26 27 28 29 30
20.4521148 26.9273190 27.9883566 15.7925567 11.9191426 5.0674371
31 32 33 34 35 36
4.7207208 -1.0048214 -8.2381733 -7.6433053 -10.2264181 -23.8004225
37 38 39 40 41 42
-12.4816686 -0.8958088 -8.5667998 -14.5701798 -15.5128079 -14.0965420
43 44 45 46 47 48
-17.6252404 -23.5605371 -21.8945413 -29.1843582 -28.0327706 -27.0525035
49 50 51 52 53 54
-18.0746007 -11.8878091 -20.5535811 -13.8988062 -8.1392908 -2.8169050
55 56 57 58 59 60
9.3723216 14.2534896 18.3494324 24.8619767 24.7770311 28.2759064
61
38.6172139
> postscript(file="/var/www/html/rcomp/tmp/6sx6q1260820143.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 -34.6687620 NA
1 -22.3828712 -34.6687620
2 -12.3123897 -22.3828712
3 -0.3231012 -12.3123897
4 0.3635345 -0.3231012
5 2.1374270 0.3635345
6 -0.3797150 2.1374270
7 -1.8808396 -0.3797150
8 0.6236152 -1.8808396
9 -5.2817965 0.6236152
10 -3.2639773 -5.2817965
11 0.4022172 -3.2639773
12 6.1557025 0.4022172
13 8.2391701 6.1557025
14 13.4444139 8.2391701
15 12.9995305 13.4444139
16 11.3694217 12.9995305
17 9.7085829 11.3694217
18 3.9119130 9.7085829
19 12.1927084 3.9119130
20 11.1596671 12.1927084
21 17.2474833 11.1596671
22 16.7461349 17.2474833
23 22.1748024 16.7461349
24 20.4521148 22.1748024
25 26.9273190 20.4521148
26 27.9883566 26.9273190
27 15.7925567 27.9883566
28 11.9191426 15.7925567
29 5.0674371 11.9191426
30 4.7207208 5.0674371
31 -1.0048214 4.7207208
32 -8.2381733 -1.0048214
33 -7.6433053 -8.2381733
34 -10.2264181 -7.6433053
35 -23.8004225 -10.2264181
36 -12.4816686 -23.8004225
37 -0.8958088 -12.4816686
38 -8.5667998 -0.8958088
39 -14.5701798 -8.5667998
40 -15.5128079 -14.5701798
41 -14.0965420 -15.5128079
42 -17.6252404 -14.0965420
43 -23.5605371 -17.6252404
44 -21.8945413 -23.5605371
45 -29.1843582 -21.8945413
46 -28.0327706 -29.1843582
47 -27.0525035 -28.0327706
48 -18.0746007 -27.0525035
49 -11.8878091 -18.0746007
50 -20.5535811 -11.8878091
51 -13.8988062 -20.5535811
52 -8.1392908 -13.8988062
53 -2.8169050 -8.1392908
54 9.3723216 -2.8169050
55 14.2534896 9.3723216
56 18.3494324 14.2534896
57 24.8619767 18.3494324
58 24.7770311 24.8619767
59 28.2759064 24.7770311
60 38.6172139 28.2759064
61 NA 38.6172139
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -22.3828712 -34.6687620
[2,] -12.3123897 -22.3828712
[3,] -0.3231012 -12.3123897
[4,] 0.3635345 -0.3231012
[5,] 2.1374270 0.3635345
[6,] -0.3797150 2.1374270
[7,] -1.8808396 -0.3797150
[8,] 0.6236152 -1.8808396
[9,] -5.2817965 0.6236152
[10,] -3.2639773 -5.2817965
[11,] 0.4022172 -3.2639773
[12,] 6.1557025 0.4022172
[13,] 8.2391701 6.1557025
[14,] 13.4444139 8.2391701
[15,] 12.9995305 13.4444139
[16,] 11.3694217 12.9995305
[17,] 9.7085829 11.3694217
[18,] 3.9119130 9.7085829
[19,] 12.1927084 3.9119130
[20,] 11.1596671 12.1927084
[21,] 17.2474833 11.1596671
[22,] 16.7461349 17.2474833
[23,] 22.1748024 16.7461349
[24,] 20.4521148 22.1748024
[25,] 26.9273190 20.4521148
[26,] 27.9883566 26.9273190
[27,] 15.7925567 27.9883566
[28,] 11.9191426 15.7925567
[29,] 5.0674371 11.9191426
[30,] 4.7207208 5.0674371
[31,] -1.0048214 4.7207208
[32,] -8.2381733 -1.0048214
[33,] -7.6433053 -8.2381733
[34,] -10.2264181 -7.6433053
[35,] -23.8004225 -10.2264181
[36,] -12.4816686 -23.8004225
[37,] -0.8958088 -12.4816686
[38,] -8.5667998 -0.8958088
[39,] -14.5701798 -8.5667998
[40,] -15.5128079 -14.5701798
[41,] -14.0965420 -15.5128079
[42,] -17.6252404 -14.0965420
[43,] -23.5605371 -17.6252404
[44,] -21.8945413 -23.5605371
[45,] -29.1843582 -21.8945413
[46,] -28.0327706 -29.1843582
[47,] -27.0525035 -28.0327706
[48,] -18.0746007 -27.0525035
[49,] -11.8878091 -18.0746007
[50,] -20.5535811 -11.8878091
[51,] -13.8988062 -20.5535811
[52,] -8.1392908 -13.8988062
[53,] -2.8169050 -8.1392908
[54,] 9.3723216 -2.8169050
[55,] 14.2534896 9.3723216
[56,] 18.3494324 14.2534896
[57,] 24.8619767 18.3494324
[58,] 24.7770311 24.8619767
[59,] 28.2759064 24.7770311
[60,] 38.6172139 28.2759064
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -22.3828712 -34.6687620
2 -12.3123897 -22.3828712
3 -0.3231012 -12.3123897
4 0.3635345 -0.3231012
5 2.1374270 0.3635345
6 -0.3797150 2.1374270
7 -1.8808396 -0.3797150
8 0.6236152 -1.8808396
9 -5.2817965 0.6236152
10 -3.2639773 -5.2817965
11 0.4022172 -3.2639773
12 6.1557025 0.4022172
13 8.2391701 6.1557025
14 13.4444139 8.2391701
15 12.9995305 13.4444139
16 11.3694217 12.9995305
17 9.7085829 11.3694217
18 3.9119130 9.7085829
19 12.1927084 3.9119130
20 11.1596671 12.1927084
21 17.2474833 11.1596671
22 16.7461349 17.2474833
23 22.1748024 16.7461349
24 20.4521148 22.1748024
25 26.9273190 20.4521148
26 27.9883566 26.9273190
27 15.7925567 27.9883566
28 11.9191426 15.7925567
29 5.0674371 11.9191426
30 4.7207208 5.0674371
31 -1.0048214 4.7207208
32 -8.2381733 -1.0048214
33 -7.6433053 -8.2381733
34 -10.2264181 -7.6433053
35 -23.8004225 -10.2264181
36 -12.4816686 -23.8004225
37 -0.8958088 -12.4816686
38 -8.5667998 -0.8958088
39 -14.5701798 -8.5667998
40 -15.5128079 -14.5701798
41 -14.0965420 -15.5128079
42 -17.6252404 -14.0965420
43 -23.5605371 -17.6252404
44 -21.8945413 -23.5605371
45 -29.1843582 -21.8945413
46 -28.0327706 -29.1843582
47 -27.0525035 -28.0327706
48 -18.0746007 -27.0525035
49 -11.8878091 -18.0746007
50 -20.5535811 -11.8878091
51 -13.8988062 -20.5535811
52 -8.1392908 -13.8988062
53 -2.8169050 -8.1392908
54 9.3723216 -2.8169050
55 14.2534896 9.3723216
56 18.3494324 14.2534896
57 24.8619767 18.3494324
58 24.7770311 24.8619767
59 28.2759064 24.7770311
60 38.6172139 28.2759064
> 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/7l9tt1260820143.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/8mvzd1260820143.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/9emhk1260820143.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/10wr6n1260820143.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/11lozl1260820143.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/122wm11260820143.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/13pwhg1260820143.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/14kw181260820143.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/15ogo21260820144.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/166t3j1260820144.tab")
+ }
>
> try(system("convert tmp/1p52h1260820143.ps tmp/1p52h1260820143.png",intern=TRUE))
character(0)
> try(system("convert tmp/212m61260820143.ps tmp/212m61260820143.png",intern=TRUE))
character(0)
> try(system("convert tmp/38oln1260820143.ps tmp/38oln1260820143.png",intern=TRUE))
character(0)
> try(system("convert tmp/4qqab1260820143.ps tmp/4qqab1260820143.png",intern=TRUE))
character(0)
> try(system("convert tmp/55yda1260820143.ps tmp/55yda1260820143.png",intern=TRUE))
character(0)
> try(system("convert tmp/6sx6q1260820143.ps tmp/6sx6q1260820143.png",intern=TRUE))
character(0)
> try(system("convert tmp/7l9tt1260820143.ps tmp/7l9tt1260820143.png",intern=TRUE))
character(0)
> try(system("convert tmp/8mvzd1260820143.ps tmp/8mvzd1260820143.png",intern=TRUE))
character(0)
> try(system("convert tmp/9emhk1260820143.ps tmp/9emhk1260820143.png",intern=TRUE))
character(0)
> try(system("convert tmp/10wr6n1260820143.ps tmp/10wr6n1260820143.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.335 1.589 3.066