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(449,0,452,0,462,0,455,0,461,0,461,0,463,0,462,0,456,0,455,0,456,0,472,0,472,0,471,0,465,0,459,0,465,0,468,0,467,0,463,0,460,0,462,0,461,0,476,0,476,0,471,0,453,0,443,0,442,0,444,0,438,0,427,0,424,0,416,0,406,0,431,0,434,0,418,0,412,0,404,0,409,0,412,1,406,1,398,1,397,1,385,1,390,1,413,1,413,1,401,1,397,1,397,1,409,1,419,1,424,1,428,1,430,1,424,1,433,1,456,1,459,1),dim=c(2,61),dimnames=list(c('Y','X'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Y','X'),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 = '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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 449 0 1 0 0 0 0 0 0 0 0 0 0
2 452 0 0 1 0 0 0 0 0 0 0 0 0
3 462 0 0 0 1 0 0 0 0 0 0 0 0
4 455 0 0 0 0 1 0 0 0 0 0 0 0
5 461 0 0 0 0 0 1 0 0 0 0 0 0
6 461 0 0 0 0 0 0 1 0 0 0 0 0
7 463 0 0 0 0 0 0 0 1 0 0 0 0
8 462 0 0 0 0 0 0 0 0 1 0 0 0
9 456 0 0 0 0 0 0 0 0 0 1 0 0
10 455 0 0 0 0 0 0 0 0 0 0 1 0
11 456 0 0 0 0 0 0 0 0 0 0 0 1
12 472 0 0 0 0 0 0 0 0 0 0 0 0
13 472 0 1 0 0 0 0 0 0 0 0 0 0
14 471 0 0 1 0 0 0 0 0 0 0 0 0
15 465 0 0 0 1 0 0 0 0 0 0 0 0
16 459 0 0 0 0 1 0 0 0 0 0 0 0
17 465 0 0 0 0 0 1 0 0 0 0 0 0
18 468 0 0 0 0 0 0 1 0 0 0 0 0
19 467 0 0 0 0 0 0 0 1 0 0 0 0
20 463 0 0 0 0 0 0 0 0 1 0 0 0
21 460 0 0 0 0 0 0 0 0 0 1 0 0
22 462 0 0 0 0 0 0 0 0 0 0 1 0
23 461 0 0 0 0 0 0 0 0 0 0 0 1
24 476 0 0 0 0 0 0 0 0 0 0 0 0
25 476 0 1 0 0 0 0 0 0 0 0 0 0
26 471 0 0 1 0 0 0 0 0 0 0 0 0
27 453 0 0 0 1 0 0 0 0 0 0 0 0
28 443 0 0 0 0 1 0 0 0 0 0 0 0
29 442 0 0 0 0 0 1 0 0 0 0 0 0
30 444 0 0 0 0 0 0 1 0 0 0 0 0
31 438 0 0 0 0 0 0 0 1 0 0 0 0
32 427 0 0 0 0 0 0 0 0 1 0 0 0
33 424 0 0 0 0 0 0 0 0 0 1 0 0
34 416 0 0 0 0 0 0 0 0 0 0 1 0
35 406 0 0 0 0 0 0 0 0 0 0 0 1
36 431 0 0 0 0 0 0 0 0 0 0 0 0
37 434 0 1 0 0 0 0 0 0 0 0 0 0
38 418 0 0 1 0 0 0 0 0 0 0 0 0
39 412 0 0 0 1 0 0 0 0 0 0 0 0
40 404 0 0 0 0 1 0 0 0 0 0 0 0
41 409 0 0 0 0 0 1 0 0 0 0 0 0
42 412 1 0 0 0 0 0 1 0 0 0 0 0
43 406 1 0 0 0 0 0 0 1 0 0 0 0
44 398 1 0 0 0 0 0 0 0 1 0 0 0
45 397 1 0 0 0 0 0 0 0 0 1 0 0
46 385 1 0 0 0 0 0 0 0 0 0 1 0
47 390 1 0 0 0 0 0 0 0 0 0 0 1
48 413 1 0 0 0 0 0 0 0 0 0 0 0
49 413 1 1 0 0 0 0 0 0 0 0 0 0
50 401 1 0 1 0 0 0 0 0 0 0 0 0
51 397 1 0 0 1 0 0 0 0 0 0 0 0
52 397 1 0 0 0 1 0 0 0 0 0 0 0
53 409 1 0 0 0 0 1 0 0 0 0 0 0
54 419 1 0 0 0 0 0 1 0 0 0 0 0
55 424 1 0 0 0 0 0 0 1 0 0 0 0
56 428 1 0 0 0 0 0 0 0 1 0 0 0
57 430 1 0 0 0 0 0 0 0 0 1 0 0
58 424 1 0 0 0 0 0 0 0 0 0 1 0
59 433 1 0 0 0 0 0 0 0 0 0 0 1
60 456 1 0 0 0 0 0 0 0 0 0 0 0
61 459 1 1 0 0 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) X M1 M2 M3 M4
464.210 -36.526 -1.535 -14.305 -19.105 -25.305
M5 M6 M7 M8 M9 M10
-19.705 -8.800 -10.000 -14.000 -16.200 -21.200
M11
-20.400
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-37.810 -14.485 6.315 14.315 32.851
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 464.210 9.704 47.836 < 2e-16 ***
X -36.526 5.855 -6.238 1.08e-07 ***
M1 -1.535 12.757 -0.120 0.9047
M2 -14.305 13.369 -1.070 0.2900
M3 -19.105 13.369 -1.429 0.1595
M4 -25.305 13.369 -1.893 0.0644 .
M5 -19.705 13.369 -1.474 0.1470
M6 -8.800 13.318 -0.661 0.5119
M7 -10.000 13.318 -0.751 0.4564
M8 -14.000 13.318 -1.051 0.2984
M9 -16.200 13.318 -1.216 0.2298
M10 -21.200 13.318 -1.592 0.1180
M11 -20.400 13.318 -1.532 0.1321
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.06 on 48 degrees of freedom
Multiple R-squared: 0.4871, Adjusted R-squared: 0.3588
F-statistic: 3.798 on 12 and 48 DF, p-value: 0.0004474
> 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,] 0.1711196892 0.3422393785 0.8288803
[2,] 0.0809937873 0.1619875746 0.9190062
[3,] 0.0372396525 0.0744793049 0.9627603
[4,] 0.0155662925 0.0311325850 0.9844337
[5,] 0.0061049785 0.0122099570 0.9938950
[6,] 0.0024069351 0.0048138702 0.9975931
[7,] 0.0012611805 0.0025223610 0.9987388
[8,] 0.0006346859 0.0012693719 0.9993653
[9,] 0.0002795451 0.0005590902 0.9997205
[10,] 0.0004266399 0.0008532797 0.9995734
[11,] 0.0006733453 0.0013466906 0.9993267
[12,] 0.0009322064 0.0018644128 0.9990678
[13,] 0.0017245622 0.0034491244 0.9982754
[14,] 0.0050325346 0.0100650692 0.9949675
[15,] 0.0087934075 0.0175868150 0.9912066
[16,] 0.0207874743 0.0415749487 0.9792125
[17,] 0.0574138796 0.1148277592 0.9425861
[18,] 0.0918919017 0.1837838035 0.9081081
[19,] 0.1653760831 0.3307521662 0.8346239
[20,] 0.2940104249 0.5880208498 0.7059896
[21,] 0.3325765114 0.6651530228 0.6674235
[22,] 0.3138184051 0.6276368103 0.6861816
[23,] 0.3433571461 0.6867142922 0.6566429
[24,] 0.3555866113 0.7111732226 0.6444134
[25,] 0.3429288475 0.6858576950 0.6570712
[26,] 0.3072759191 0.6145518382 0.6927241
[27,] 0.2100596021 0.4201192042 0.7899404
[28,] 0.1409363559 0.2818727119 0.8590636
[29,] 0.1031873982 0.2063747964 0.8968126
[30,] 0.0746860518 0.1493721037 0.9253139
> postscript(file="/var/www/html/rcomp/tmp/1jjyp1258658318.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/20ipl1258658318.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/3zpon1258658318.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/4map91258658318.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/53kh01258658318.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
-13.6752577 2.0948454 16.8948454 16.0948454 16.4948454 5.5896907
7 8 9 10 11 12
8.7896907 11.7896907 7.9896907 11.9896907 12.1896907 7.7896907
13 14 15 16 17 18
9.3247423 21.0948454 19.8948454 20.0948454 20.4948454 12.5896907
19 20 21 22 23 24
12.7896907 12.7896907 11.9896907 18.9896907 17.1896907 11.7896907
25 26 27 28 29 30
13.3247423 21.0948454 7.8948454 4.0948454 -2.5051546 -11.4103093
31 32 33 34 35 36
-16.2103093 -23.2103093 -24.0103093 -27.0103093 -37.8103093 -33.2103093
37 38 39 40 41 42
-28.6752577 -31.9051546 -33.1051546 -34.9051546 -35.5051546 -6.8845361
43 44 45 46 47 48
-11.6845361 -15.6845361 -14.4845361 -21.4845361 -17.2845361 -14.6845361
49 50 51 52 53 54
-13.1494845 -12.3793814 -11.5793814 -5.3793814 1.0206186 0.1154639
55 56 57 58 59 60
6.3154639 14.3154639 18.5154639 17.5154639 25.7154639 28.3154639
61
32.8505155
> postscript(file="/var/www/html/rcomp/tmp/6gnz01258658318.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 -13.6752577 NA
1 2.0948454 -13.6752577
2 16.8948454 2.0948454
3 16.0948454 16.8948454
4 16.4948454 16.0948454
5 5.5896907 16.4948454
6 8.7896907 5.5896907
7 11.7896907 8.7896907
8 7.9896907 11.7896907
9 11.9896907 7.9896907
10 12.1896907 11.9896907
11 7.7896907 12.1896907
12 9.3247423 7.7896907
13 21.0948454 9.3247423
14 19.8948454 21.0948454
15 20.0948454 19.8948454
16 20.4948454 20.0948454
17 12.5896907 20.4948454
18 12.7896907 12.5896907
19 12.7896907 12.7896907
20 11.9896907 12.7896907
21 18.9896907 11.9896907
22 17.1896907 18.9896907
23 11.7896907 17.1896907
24 13.3247423 11.7896907
25 21.0948454 13.3247423
26 7.8948454 21.0948454
27 4.0948454 7.8948454
28 -2.5051546 4.0948454
29 -11.4103093 -2.5051546
30 -16.2103093 -11.4103093
31 -23.2103093 -16.2103093
32 -24.0103093 -23.2103093
33 -27.0103093 -24.0103093
34 -37.8103093 -27.0103093
35 -33.2103093 -37.8103093
36 -28.6752577 -33.2103093
37 -31.9051546 -28.6752577
38 -33.1051546 -31.9051546
39 -34.9051546 -33.1051546
40 -35.5051546 -34.9051546
41 -6.8845361 -35.5051546
42 -11.6845361 -6.8845361
43 -15.6845361 -11.6845361
44 -14.4845361 -15.6845361
45 -21.4845361 -14.4845361
46 -17.2845361 -21.4845361
47 -14.6845361 -17.2845361
48 -13.1494845 -14.6845361
49 -12.3793814 -13.1494845
50 -11.5793814 -12.3793814
51 -5.3793814 -11.5793814
52 1.0206186 -5.3793814
53 0.1154639 1.0206186
54 6.3154639 0.1154639
55 14.3154639 6.3154639
56 18.5154639 14.3154639
57 17.5154639 18.5154639
58 25.7154639 17.5154639
59 28.3154639 25.7154639
60 32.8505155 28.3154639
61 NA 32.8505155
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 2.0948454 -13.6752577
[2,] 16.8948454 2.0948454
[3,] 16.0948454 16.8948454
[4,] 16.4948454 16.0948454
[5,] 5.5896907 16.4948454
[6,] 8.7896907 5.5896907
[7,] 11.7896907 8.7896907
[8,] 7.9896907 11.7896907
[9,] 11.9896907 7.9896907
[10,] 12.1896907 11.9896907
[11,] 7.7896907 12.1896907
[12,] 9.3247423 7.7896907
[13,] 21.0948454 9.3247423
[14,] 19.8948454 21.0948454
[15,] 20.0948454 19.8948454
[16,] 20.4948454 20.0948454
[17,] 12.5896907 20.4948454
[18,] 12.7896907 12.5896907
[19,] 12.7896907 12.7896907
[20,] 11.9896907 12.7896907
[21,] 18.9896907 11.9896907
[22,] 17.1896907 18.9896907
[23,] 11.7896907 17.1896907
[24,] 13.3247423 11.7896907
[25,] 21.0948454 13.3247423
[26,] 7.8948454 21.0948454
[27,] 4.0948454 7.8948454
[28,] -2.5051546 4.0948454
[29,] -11.4103093 -2.5051546
[30,] -16.2103093 -11.4103093
[31,] -23.2103093 -16.2103093
[32,] -24.0103093 -23.2103093
[33,] -27.0103093 -24.0103093
[34,] -37.8103093 -27.0103093
[35,] -33.2103093 -37.8103093
[36,] -28.6752577 -33.2103093
[37,] -31.9051546 -28.6752577
[38,] -33.1051546 -31.9051546
[39,] -34.9051546 -33.1051546
[40,] -35.5051546 -34.9051546
[41,] -6.8845361 -35.5051546
[42,] -11.6845361 -6.8845361
[43,] -15.6845361 -11.6845361
[44,] -14.4845361 -15.6845361
[45,] -21.4845361 -14.4845361
[46,] -17.2845361 -21.4845361
[47,] -14.6845361 -17.2845361
[48,] -13.1494845 -14.6845361
[49,] -12.3793814 -13.1494845
[50,] -11.5793814 -12.3793814
[51,] -5.3793814 -11.5793814
[52,] 1.0206186 -5.3793814
[53,] 0.1154639 1.0206186
[54,] 6.3154639 0.1154639
[55,] 14.3154639 6.3154639
[56,] 18.5154639 14.3154639
[57,] 17.5154639 18.5154639
[58,] 25.7154639 17.5154639
[59,] 28.3154639 25.7154639
[60,] 32.8505155 28.3154639
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 2.0948454 -13.6752577
2 16.8948454 2.0948454
3 16.0948454 16.8948454
4 16.4948454 16.0948454
5 5.5896907 16.4948454
6 8.7896907 5.5896907
7 11.7896907 8.7896907
8 7.9896907 11.7896907
9 11.9896907 7.9896907
10 12.1896907 11.9896907
11 7.7896907 12.1896907
12 9.3247423 7.7896907
13 21.0948454 9.3247423
14 19.8948454 21.0948454
15 20.0948454 19.8948454
16 20.4948454 20.0948454
17 12.5896907 20.4948454
18 12.7896907 12.5896907
19 12.7896907 12.7896907
20 11.9896907 12.7896907
21 18.9896907 11.9896907
22 17.1896907 18.9896907
23 11.7896907 17.1896907
24 13.3247423 11.7896907
25 21.0948454 13.3247423
26 7.8948454 21.0948454
27 4.0948454 7.8948454
28 -2.5051546 4.0948454
29 -11.4103093 -2.5051546
30 -16.2103093 -11.4103093
31 -23.2103093 -16.2103093
32 -24.0103093 -23.2103093
33 -27.0103093 -24.0103093
34 -37.8103093 -27.0103093
35 -33.2103093 -37.8103093
36 -28.6752577 -33.2103093
37 -31.9051546 -28.6752577
38 -33.1051546 -31.9051546
39 -34.9051546 -33.1051546
40 -35.5051546 -34.9051546
41 -6.8845361 -35.5051546
42 -11.6845361 -6.8845361
43 -15.6845361 -11.6845361
44 -14.4845361 -15.6845361
45 -21.4845361 -14.4845361
46 -17.2845361 -21.4845361
47 -14.6845361 -17.2845361
48 -13.1494845 -14.6845361
49 -12.3793814 -13.1494845
50 -11.5793814 -12.3793814
51 -5.3793814 -11.5793814
52 1.0206186 -5.3793814
53 0.1154639 1.0206186
54 6.3154639 0.1154639
55 14.3154639 6.3154639
56 18.5154639 14.3154639
57 17.5154639 18.5154639
58 25.7154639 17.5154639
59 28.3154639 25.7154639
60 32.8505155 28.3154639
> 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/7mo5w1258658318.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/8d7pw1258658318.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/9opz81258658318.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/10rwfb1258658318.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/11x1zt1258658318.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/12ifd91258658318.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/13kyfh1258658318.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/14o7uu1258658318.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/155ajg1258658318.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/16pn1d1258658318.tab")
+ }
>
> system("convert tmp/1jjyp1258658318.ps tmp/1jjyp1258658318.png")
> system("convert tmp/20ipl1258658318.ps tmp/20ipl1258658318.png")
> system("convert tmp/3zpon1258658318.ps tmp/3zpon1258658318.png")
> system("convert tmp/4map91258658318.ps tmp/4map91258658318.png")
> system("convert tmp/53kh01258658318.ps tmp/53kh01258658318.png")
> system("convert tmp/6gnz01258658318.ps tmp/6gnz01258658318.png")
> system("convert tmp/7mo5w1258658318.ps tmp/7mo5w1258658318.png")
> system("convert tmp/8d7pw1258658318.ps tmp/8d7pw1258658318.png")
> system("convert tmp/9opz81258658318.ps tmp/9opz81258658318.png")
> system("convert tmp/10rwfb1258658318.ps tmp/10rwfb1258658318.png")
>
>
> proc.time()
user system elapsed
2.436 1.582 2.826