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(462,1919,455,1911,461,1870,461,2263,463,1802,462,1863,456,1989,455,2197,456,2409,472,2502,472,2593,471,2598,465,2053,459,2213,465,2238,468,2359,467,2151,463,2474,460,3079,462,2312,461,2565,476,1972,476,2484,471,2202,453,2151,443,1976,442,2012,444,2114,438,1772,427,1957,424,2070,416,1990,406,2182,431,2008,434,1916,418,2397,412,2114,404,1778,409,1641,412,2186,406,1773,398,1785,397,2217,385,2153,390,1895,413,2475,413,1793,401,2308,397,2051,397,1898,409,2142,419,1874,424,1560,428,1808,430,1575,424,1525,433,1997,456,1753,459,1623,446,2251,441,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 462 1919 1 0 0 0 0 0 0 0 0 0 0 1
2 455 1911 0 1 0 0 0 0 0 0 0 0 0 2
3 461 1870 0 0 1 0 0 0 0 0 0 0 0 3
4 461 2263 0 0 0 1 0 0 0 0 0 0 0 4
5 463 1802 0 0 0 0 1 0 0 0 0 0 0 5
6 462 1863 0 0 0 0 0 1 0 0 0 0 0 6
7 456 1989 0 0 0 0 0 0 1 0 0 0 0 7
8 455 2197 0 0 0 0 0 0 0 1 0 0 0 8
9 456 2409 0 0 0 0 0 0 0 0 1 0 0 9
10 472 2502 0 0 0 0 0 0 0 0 0 1 0 10
11 472 2593 0 0 0 0 0 0 0 0 0 0 1 11
12 471 2598 0 0 0 0 0 0 0 0 0 0 0 12
13 465 2053 1 0 0 0 0 0 0 0 0 0 0 13
14 459 2213 0 1 0 0 0 0 0 0 0 0 0 14
15 465 2238 0 0 1 0 0 0 0 0 0 0 0 15
16 468 2359 0 0 0 1 0 0 0 0 0 0 0 16
17 467 2151 0 0 0 0 1 0 0 0 0 0 0 17
18 463 2474 0 0 0 0 0 1 0 0 0 0 0 18
19 460 3079 0 0 0 0 0 0 1 0 0 0 0 19
20 462 2312 0 0 0 0 0 0 0 1 0 0 0 20
21 461 2565 0 0 0 0 0 0 0 0 1 0 0 21
22 476 1972 0 0 0 0 0 0 0 0 0 1 0 22
23 476 2484 0 0 0 0 0 0 0 0 0 0 1 23
24 471 2202 0 0 0 0 0 0 0 0 0 0 0 24
25 453 2151 1 0 0 0 0 0 0 0 0 0 0 25
26 443 1976 0 1 0 0 0 0 0 0 0 0 0 26
27 442 2012 0 0 1 0 0 0 0 0 0 0 0 27
28 444 2114 0 0 0 1 0 0 0 0 0 0 0 28
29 438 1772 0 0 0 0 1 0 0 0 0 0 0 29
30 427 1957 0 0 0 0 0 1 0 0 0 0 0 30
31 424 2070 0 0 0 0 0 0 1 0 0 0 0 31
32 416 1990 0 0 0 0 0 0 0 1 0 0 0 32
33 406 2182 0 0 0 0 0 0 0 0 1 0 0 33
34 431 2008 0 0 0 0 0 0 0 0 0 1 0 34
35 434 1916 0 0 0 0 0 0 0 0 0 0 1 35
36 418 2397 0 0 0 0 0 0 0 0 0 0 0 36
37 412 2114 1 0 0 0 0 0 0 0 0 0 0 37
38 404 1778 0 1 0 0 0 0 0 0 0 0 0 38
39 409 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 406 1773 0 0 0 0 1 0 0 0 0 0 0 41
42 398 1785 0 0 0 0 0 1 0 0 0 0 0 42
43 397 2217 0 0 0 0 0 0 1 0 0 0 0 43
44 385 2153 0 0 0 0 0 0 0 1 0 0 0 44
45 390 1895 0 0 0 0 0 0 0 0 1 0 0 45
46 413 2475 0 0 0 0 0 0 0 0 0 1 0 46
47 413 1793 0 0 0 0 0 0 0 0 0 0 1 47
48 401 2308 0 0 0 0 0 0 0 0 0 0 0 48
49 397 2051 1 0 0 0 0 0 0 0 0 0 0 49
50 397 1898 0 1 0 0 0 0 0 0 0 0 0 50
51 409 2142 0 0 1 0 0 0 0 0 0 0 0 51
52 419 1874 0 0 0 1 0 0 0 0 0 0 0 52
53 424 1560 0 0 0 0 1 0 0 0 0 0 0 53
54 428 1808 0 0 0 0 0 1 0 0 0 0 0 54
55 430 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 433 1997 0 0 0 0 0 0 0 0 1 0 0 57
58 456 1753 0 0 0 0 0 0 0 0 0 1 0 58
59 459 1623 0 0 0 0 0 0 0 0 0 0 1 59
60 446 2251 0 0 0 0 0 0 0 0 0 0 0 60
61 441 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
4.780e+02 3.370e-04 -8.151e+00 -2.005e+01 -1.342e+01 -8.843e+00
M5 M6 M7 M8 M9 M10
-8.888e+00 -1.191e+01 -1.314e+01 -1.705e+01 -1.527e+01 6.193e+00
M11 t
8.452e+00 -1.039e+00
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-30.978 -15.623 -1.087 14.671 33.869
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 477.993761 31.666012 15.095 < 2e-16 ***
bvg 0.000337 0.011353 0.030 0.976
M1 -8.150826 12.732764 -0.640 0.525
M2 -20.051592 13.754834 -1.458 0.152
M3 -13.421647 13.595873 -0.987 0.329
M4 -8.843337 12.950613 -0.683 0.498
M5 -8.887673 14.304897 -0.621 0.537
M6 -11.905050 13.475594 -0.883 0.381
M7 -13.136851 12.801573 -1.026 0.310
M8 -17.047586 13.184982 -1.293 0.202
M9 -15.267793 12.714327 -1.201 0.236
M10 6.193498 12.827693 0.483 0.631
M11 8.452294 12.956477 0.652 0.517
t -1.038506 0.172032 -6.037 2.37e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 19.87 on 47 degrees of freedom
Multiple R-squared: 0.5502, Adjusted R-squared: 0.4257
F-statistic: 4.422 on 13 and 47 DF, p-value: 7.768e-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.405426e-04 2.810852e-04 0.9998594574
[2,] 6.481734e-06 1.296347e-05 0.9999935183
[3,] 1.614455e-06 3.228910e-06 0.9999983855
[4,] 2.202815e-07 4.405629e-07 0.9999997797
[5,] 1.738755e-08 3.477511e-08 0.9999999826
[6,] 1.797744e-09 3.595488e-09 0.9999999982
[7,] 2.336947e-10 4.673894e-10 0.9999999998
[8,] 2.562392e-10 5.124783e-10 0.9999999997
[9,] 1.281503e-06 2.563006e-06 0.9999987185
[10,] 1.746068e-05 3.492137e-05 0.9999825393
[11,] 2.098107e-04 4.196213e-04 0.9997901893
[12,] 4.801066e-04 9.602132e-04 0.9995198934
[13,] 1.687925e-03 3.375850e-03 0.9983120748
[14,] 7.889516e-03 1.577903e-02 0.9921104842
[15,] 1.188910e-02 2.377819e-02 0.9881109034
[16,] 4.274069e-02 8.548137e-02 0.9572593135
[17,] 1.074952e-01 2.149904e-01 0.8925047834
[18,] 1.104116e-01 2.208231e-01 0.8895884414
[19,] 1.264628e-01 2.529255e-01 0.8735372343
[20,] 2.809443e-01 5.618887e-01 0.7190556695
[21,] 4.590694e-01 9.181387e-01 0.5409306264
[22,] 6.076020e-01 7.847960e-01 0.3923979773
[23,] 7.260636e-01 5.478728e-01 0.2739364100
[24,] 9.295095e-01 1.409811e-01 0.0704905350
[25,] 9.900956e-01 1.980887e-02 0.0099044365
[26,] 9.981499e-01 3.700100e-03 0.0018500502
[27,] 9.993595e-01 1.281086e-03 0.0006405431
[28,] 9.991609e-01 1.678158e-03 0.0008390791
> postscript(file="/var/www/html/rcomp/tmp/1fnum1258742311.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/2kctz1258742311.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/3y0ue1258742311.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/4taf91258742311.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/5he641258742311.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 7
-7.451223 -1.509254 -1.086874 -4.759137 -1.520916 1.514407 -2.257753
8 9 10 11 12 13 14
1.621382 1.808641 -2.645488 -3.896450 4.592665 7.965689 14.851033
15 16 17 18 19 20 21
15.251168 14.670582 14.823530 14.770547 13.836941 21.044697 19.218138
22 23 24 25 26 27 28
13.995222 12.602364 17.188212 8.394734 11.392989 4.789417 3.215235
29 30 31 32 33 34 35
-1.586653 -8.593124 -9.360902 -12.384698 -23.190697 -18.554836 -16.744118
36 37 38 39 40 41 42
-23.415437 -20.130720 -15.078200 -15.623463 -16.346957 -21.124914 -25.073076
43 44 45 46 47 48 49
-23.948373 -30.977561 -26.631889 -24.250161 -25.240585 -27.923364 -22.647410
50 51 52 53 54 55 56
-9.656570 -3.330248 3.220277 9.408952 17.381247 21.730087 20.696180
57 58 59 60 61
28.795808 31.455263 33.278789 29.557923 33.868930
> postscript(file="/var/www/html/rcomp/tmp/6q0201258742311.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 -7.451223 NA
1 -1.509254 -7.451223
2 -1.086874 -1.509254
3 -4.759137 -1.086874
4 -1.520916 -4.759137
5 1.514407 -1.520916
6 -2.257753 1.514407
7 1.621382 -2.257753
8 1.808641 1.621382
9 -2.645488 1.808641
10 -3.896450 -2.645488
11 4.592665 -3.896450
12 7.965689 4.592665
13 14.851033 7.965689
14 15.251168 14.851033
15 14.670582 15.251168
16 14.823530 14.670582
17 14.770547 14.823530
18 13.836941 14.770547
19 21.044697 13.836941
20 19.218138 21.044697
21 13.995222 19.218138
22 12.602364 13.995222
23 17.188212 12.602364
24 8.394734 17.188212
25 11.392989 8.394734
26 4.789417 11.392989
27 3.215235 4.789417
28 -1.586653 3.215235
29 -8.593124 -1.586653
30 -9.360902 -8.593124
31 -12.384698 -9.360902
32 -23.190697 -12.384698
33 -18.554836 -23.190697
34 -16.744118 -18.554836
35 -23.415437 -16.744118
36 -20.130720 -23.415437
37 -15.078200 -20.130720
38 -15.623463 -15.078200
39 -16.346957 -15.623463
40 -21.124914 -16.346957
41 -25.073076 -21.124914
42 -23.948373 -25.073076
43 -30.977561 -23.948373
44 -26.631889 -30.977561
45 -24.250161 -26.631889
46 -25.240585 -24.250161
47 -27.923364 -25.240585
48 -22.647410 -27.923364
49 -9.656570 -22.647410
50 -3.330248 -9.656570
51 3.220277 -3.330248
52 9.408952 3.220277
53 17.381247 9.408952
54 21.730087 17.381247
55 20.696180 21.730087
56 28.795808 20.696180
57 31.455263 28.795808
58 33.278789 31.455263
59 29.557923 33.278789
60 33.868930 29.557923
61 NA 33.868930
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1.509254 -7.451223
[2,] -1.086874 -1.509254
[3,] -4.759137 -1.086874
[4,] -1.520916 -4.759137
[5,] 1.514407 -1.520916
[6,] -2.257753 1.514407
[7,] 1.621382 -2.257753
[8,] 1.808641 1.621382
[9,] -2.645488 1.808641
[10,] -3.896450 -2.645488
[11,] 4.592665 -3.896450
[12,] 7.965689 4.592665
[13,] 14.851033 7.965689
[14,] 15.251168 14.851033
[15,] 14.670582 15.251168
[16,] 14.823530 14.670582
[17,] 14.770547 14.823530
[18,] 13.836941 14.770547
[19,] 21.044697 13.836941
[20,] 19.218138 21.044697
[21,] 13.995222 19.218138
[22,] 12.602364 13.995222
[23,] 17.188212 12.602364
[24,] 8.394734 17.188212
[25,] 11.392989 8.394734
[26,] 4.789417 11.392989
[27,] 3.215235 4.789417
[28,] -1.586653 3.215235
[29,] -8.593124 -1.586653
[30,] -9.360902 -8.593124
[31,] -12.384698 -9.360902
[32,] -23.190697 -12.384698
[33,] -18.554836 -23.190697
[34,] -16.744118 -18.554836
[35,] -23.415437 -16.744118
[36,] -20.130720 -23.415437
[37,] -15.078200 -20.130720
[38,] -15.623463 -15.078200
[39,] -16.346957 -15.623463
[40,] -21.124914 -16.346957
[41,] -25.073076 -21.124914
[42,] -23.948373 -25.073076
[43,] -30.977561 -23.948373
[44,] -26.631889 -30.977561
[45,] -24.250161 -26.631889
[46,] -25.240585 -24.250161
[47,] -27.923364 -25.240585
[48,] -22.647410 -27.923364
[49,] -9.656570 -22.647410
[50,] -3.330248 -9.656570
[51,] 3.220277 -3.330248
[52,] 9.408952 3.220277
[53,] 17.381247 9.408952
[54,] 21.730087 17.381247
[55,] 20.696180 21.730087
[56,] 28.795808 20.696180
[57,] 31.455263 28.795808
[58,] 33.278789 31.455263
[59,] 29.557923 33.278789
[60,] 33.868930 29.557923
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1.509254 -7.451223
2 -1.086874 -1.509254
3 -4.759137 -1.086874
4 -1.520916 -4.759137
5 1.514407 -1.520916
6 -2.257753 1.514407
7 1.621382 -2.257753
8 1.808641 1.621382
9 -2.645488 1.808641
10 -3.896450 -2.645488
11 4.592665 -3.896450
12 7.965689 4.592665
13 14.851033 7.965689
14 15.251168 14.851033
15 14.670582 15.251168
16 14.823530 14.670582
17 14.770547 14.823530
18 13.836941 14.770547
19 21.044697 13.836941
20 19.218138 21.044697
21 13.995222 19.218138
22 12.602364 13.995222
23 17.188212 12.602364
24 8.394734 17.188212
25 11.392989 8.394734
26 4.789417 11.392989
27 3.215235 4.789417
28 -1.586653 3.215235
29 -8.593124 -1.586653
30 -9.360902 -8.593124
31 -12.384698 -9.360902
32 -23.190697 -12.384698
33 -18.554836 -23.190697
34 -16.744118 -18.554836
35 -23.415437 -16.744118
36 -20.130720 -23.415437
37 -15.078200 -20.130720
38 -15.623463 -15.078200
39 -16.346957 -15.623463
40 -21.124914 -16.346957
41 -25.073076 -21.124914
42 -23.948373 -25.073076
43 -30.977561 -23.948373
44 -26.631889 -30.977561
45 -24.250161 -26.631889
46 -25.240585 -24.250161
47 -27.923364 -25.240585
48 -22.647410 -27.923364
49 -9.656570 -22.647410
50 -3.330248 -9.656570
51 3.220277 -3.330248
52 9.408952 3.220277
53 17.381247 9.408952
54 21.730087 17.381247
55 20.696180 21.730087
56 28.795808 20.696180
57 31.455263 28.795808
58 33.278789 31.455263
59 29.557923 33.278789
60 33.868930 29.557923
> 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/72jbv1258742311.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/8bd201258742311.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/97vbm1258742311.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/10noey1258742311.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/11uezn1258742311.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/12htlp1258742311.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/13tj3x1258742311.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/14e6vi1258742311.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/157eh01258742311.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/16euux1258742311.tab")
+ }
>
> system("convert tmp/1fnum1258742311.ps tmp/1fnum1258742311.png")
> system("convert tmp/2kctz1258742311.ps tmp/2kctz1258742311.png")
> system("convert tmp/3y0ue1258742311.ps tmp/3y0ue1258742311.png")
> system("convert tmp/4taf91258742311.ps tmp/4taf91258742311.png")
> system("convert tmp/5he641258742311.ps tmp/5he641258742311.png")
> system("convert tmp/6q0201258742311.ps tmp/6q0201258742311.png")
> system("convert tmp/72jbv1258742311.ps tmp/72jbv1258742311.png")
> system("convert tmp/8bd201258742311.ps tmp/8bd201258742311.png")
> system("convert tmp/97vbm1258742311.ps tmp/97vbm1258742311.png")
> system("convert tmp/10noey1258742311.ps tmp/10noey1258742311.png")
>
>
> proc.time()
user system elapsed
2.409 1.551 2.810