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(130,0,136.7,0,138.1,0,139.5,0,140.4,0,144.6,0,151.4,0,147.9,0,141.5,0,143.8,0,143.6,0,150.5,0,150.1,0,154.9,0,162.1,0,176.7,0,186.6,0,194.8,0,196.3,0,228.8,0,267.2,0,237.2,0,254.7,0,258.2,0,257.9,0,269.6,0,266.9,0,269.6,0,253.9,0,258.6,0,274.2,0,301.5,0,304.5,0,285.1,0,287.7,0,265.5,0,264.1,0,276.1,0,258.9,0,239.1,0,250.1,1,276.8,1,297.6,1,295.4,1,283,1,275.8,1,279.7,1,254.6,1,234.6,1,176.9,1,148.1,1,122.7,1,124.9,1,121.6,1,128.4,1,144.5,1,151.8,1,167.1,1,173.8,1,203.7,1,199.8,1),dim=c(2,61),dimnames=list(c('Y(t)','X(t)'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Y(t)','X(t)'),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(t) X(t) M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 130.0 0 1 0 0 0 0 0 0 0 0 0 0
2 136.7 0 0 1 0 0 0 0 0 0 0 0 0
3 138.1 0 0 0 1 0 0 0 0 0 0 0 0
4 139.5 0 0 0 0 1 0 0 0 0 0 0 0
5 140.4 0 0 0 0 0 1 0 0 0 0 0 0
6 144.6 0 0 0 0 0 0 1 0 0 0 0 0
7 151.4 0 0 0 0 0 0 0 1 0 0 0 0
8 147.9 0 0 0 0 0 0 0 0 1 0 0 0
9 141.5 0 0 0 0 0 0 0 0 0 1 0 0
10 143.8 0 0 0 0 0 0 0 0 0 0 1 0
11 143.6 0 0 0 0 0 0 0 0 0 0 0 1
12 150.5 0 0 0 0 0 0 0 0 0 0 0 0
13 150.1 0 1 0 0 0 0 0 0 0 0 0 0
14 154.9 0 0 1 0 0 0 0 0 0 0 0 0
15 162.1 0 0 0 1 0 0 0 0 0 0 0 0
16 176.7 0 0 0 0 1 0 0 0 0 0 0 0
17 186.6 0 0 0 0 0 1 0 0 0 0 0 0
18 194.8 0 0 0 0 0 0 1 0 0 0 0 0
19 196.3 0 0 0 0 0 0 0 1 0 0 0 0
20 228.8 0 0 0 0 0 0 0 0 1 0 0 0
21 267.2 0 0 0 0 0 0 0 0 0 1 0 0
22 237.2 0 0 0 0 0 0 0 0 0 0 1 0
23 254.7 0 0 0 0 0 0 0 0 0 0 0 1
24 258.2 0 0 0 0 0 0 0 0 0 0 0 0
25 257.9 0 1 0 0 0 0 0 0 0 0 0 0
26 269.6 0 0 1 0 0 0 0 0 0 0 0 0
27 266.9 0 0 0 1 0 0 0 0 0 0 0 0
28 269.6 0 0 0 0 1 0 0 0 0 0 0 0
29 253.9 0 0 0 0 0 1 0 0 0 0 0 0
30 258.6 0 0 0 0 0 0 1 0 0 0 0 0
31 274.2 0 0 0 0 0 0 0 1 0 0 0 0
32 301.5 0 0 0 0 0 0 0 0 1 0 0 0
33 304.5 0 0 0 0 0 0 0 0 0 1 0 0
34 285.1 0 0 0 0 0 0 0 0 0 0 1 0
35 287.7 0 0 0 0 0 0 0 0 0 0 0 1
36 265.5 0 0 0 0 0 0 0 0 0 0 0 0
37 264.1 0 1 0 0 0 0 0 0 0 0 0 0
38 276.1 0 0 1 0 0 0 0 0 0 0 0 0
39 258.9 0 0 0 1 0 0 0 0 0 0 0 0
40 239.1 0 0 0 0 1 0 0 0 0 0 0 0
41 250.1 1 0 0 0 0 1 0 0 0 0 0 0
42 276.8 1 0 0 0 0 0 1 0 0 0 0 0
43 297.6 1 0 0 0 0 0 0 1 0 0 0 0
44 295.4 1 0 0 0 0 0 0 0 1 0 0 0
45 283.0 1 0 0 0 0 0 0 0 0 1 0 0
46 275.8 1 0 0 0 0 0 0 0 0 0 1 0
47 279.7 1 0 0 0 0 0 0 0 0 0 0 1
48 254.6 1 0 0 0 0 0 0 0 0 0 0 0
49 234.6 1 1 0 0 0 0 0 0 0 0 0 0
50 176.9 1 0 1 0 0 0 0 0 0 0 0 0
51 148.1 1 0 0 1 0 0 0 0 0 0 0 0
52 122.7 1 0 0 0 1 0 0 0 0 0 0 0
53 124.9 1 0 0 0 0 1 0 0 0 0 0 0
54 121.6 1 0 0 0 0 0 1 0 0 0 0 0
55 128.4 1 0 0 0 0 0 0 1 0 0 0 0
56 144.5 1 0 0 0 0 0 0 0 1 0 0 0
57 151.8 1 0 0 0 0 0 0 0 0 1 0 0
58 167.1 1 0 0 0 0 0 0 0 0 0 1 0
59 173.8 1 0 0 0 0 0 0 0 0 0 0 1
60 203.7 1 0 0 0 0 0 0 0 0 0 0 0
61 199.8 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(t)` M1 M2 M3 M4
230.92 -11.05 -21.15 -25.87 -33.89 -39.19
M5 M6 M7 M8 M9 M10
-35.32 -27.22 -16.92 -2.88 3.10 -4.70
M11
1.40
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-92.521 -58.930 0.759 58.432 94.652
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 230.92 30.31 7.619 8.32e-10 ***
`X(t)` -11.05 18.03 -0.613 0.543
M1 -21.15 39.88 -0.530 0.598
M2 -25.87 41.79 -0.619 0.539
M3 -33.89 41.79 -0.811 0.421
M4 -39.19 41.79 -0.938 0.353
M5 -35.32 41.63 -0.848 0.400
M6 -27.22 41.63 -0.654 0.516
M7 -16.92 41.63 -0.406 0.686
M8 -2.88 41.63 -0.069 0.945
M9 3.10 41.63 0.074 0.941
M10 -4.70 41.63 -0.113 0.911
M11 1.40 41.63 0.034 0.973
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 65.82 on 48 degrees of freedom
Multiple R-squared: 0.0634, Adjusted R-squared: -0.1707
F-statistic: 0.2708 on 12 and 48 DF, p-value: 0.9913
> 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.07567945 0.1513589 0.9243205
[2,] 0.06945471 0.1389094 0.9305453
[3,] 0.06693101 0.1338620 0.9330690
[4,] 0.05731525 0.1146305 0.9426847
[5,] 0.10889150 0.2177830 0.8911085
[6,] 0.29879462 0.5975892 0.7012054
[7,] 0.35241196 0.7048239 0.6475880
[8,] 0.43611303 0.8722261 0.5638870
[9,] 0.48310202 0.9662040 0.5168980
[10,] 0.57828056 0.8434389 0.4217194
[11,] 0.64908674 0.7018265 0.3509133
[12,] 0.68240156 0.6351969 0.3175984
[13,] 0.70566167 0.5886767 0.2943383
[14,] 0.67437695 0.6512461 0.3256230
[15,] 0.63737850 0.7252430 0.3626215
[16,] 0.60883080 0.7823384 0.3911692
[17,] 0.58883382 0.8223324 0.4111662
[18,] 0.54766122 0.9046776 0.4523388
[19,] 0.49517867 0.9903573 0.5048213
[20,] 0.43467822 0.8693564 0.5653218
[21,] 0.36892096 0.7378419 0.6310790
[22,] 0.32516370 0.6503274 0.6748363
[23,] 0.26676766 0.5335353 0.7332323
[24,] 0.20083856 0.4016771 0.7991614
[25,] 0.13657466 0.2731493 0.8634253
[26,] 0.12399929 0.2479986 0.8760007
[27,] 0.14547735 0.2909547 0.8545226
[28,] 0.21463212 0.4292642 0.7853679
[29,] 0.29522603 0.5904521 0.7047740
[30,] 0.37231650 0.7446330 0.6276835
> postscript(file="/var/www/html/rcomp/tmp/1rs0l1259397052.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/2wjbv1259397052.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/3f0pi1259397052.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/4j70k1259397052.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/5hmon1259397052.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 8
-79.7675 -68.3505 -58.9305 -52.2305 -55.2010 -59.1010 -62.6010 -80.1410
9 10 11 12 13 14 15 16
-92.5210 -82.4210 -88.7210 -80.4210 -59.6675 -50.1505 -34.9305 -15.0305
17 18 19 20 21 22 23 24
-9.0010 -8.9010 -17.7010 0.7590 33.1790 10.9790 22.3790 27.2790
25 26 27 28 29 30 31 32
48.1325 64.5495 69.8695 77.8695 58.2990 54.8990 60.1990 73.4590
33 34 35 36 37 38 39 40
70.4790 58.8790 55.3790 34.5790 54.3325 71.0495 61.8695 47.3695
41 42 43 44 45 46 47 48
65.5515 84.1515 94.6515 78.4115 60.0315 60.6315 58.4315 34.7315
49 50 51 52 53 54 55 56
35.8850 -17.0980 -37.8780 -57.9780 -59.6485 -71.0485 -74.5485 -72.4885
57 58 59 60 61
-71.1685 -48.0685 -47.4685 -16.1685 1.0850
> postscript(file="/var/www/html/rcomp/tmp/60c6i1259397052.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 -79.7675 NA
1 -68.3505 -79.7675
2 -58.9305 -68.3505
3 -52.2305 -58.9305
4 -55.2010 -52.2305
5 -59.1010 -55.2010
6 -62.6010 -59.1010
7 -80.1410 -62.6010
8 -92.5210 -80.1410
9 -82.4210 -92.5210
10 -88.7210 -82.4210
11 -80.4210 -88.7210
12 -59.6675 -80.4210
13 -50.1505 -59.6675
14 -34.9305 -50.1505
15 -15.0305 -34.9305
16 -9.0010 -15.0305
17 -8.9010 -9.0010
18 -17.7010 -8.9010
19 0.7590 -17.7010
20 33.1790 0.7590
21 10.9790 33.1790
22 22.3790 10.9790
23 27.2790 22.3790
24 48.1325 27.2790
25 64.5495 48.1325
26 69.8695 64.5495
27 77.8695 69.8695
28 58.2990 77.8695
29 54.8990 58.2990
30 60.1990 54.8990
31 73.4590 60.1990
32 70.4790 73.4590
33 58.8790 70.4790
34 55.3790 58.8790
35 34.5790 55.3790
36 54.3325 34.5790
37 71.0495 54.3325
38 61.8695 71.0495
39 47.3695 61.8695
40 65.5515 47.3695
41 84.1515 65.5515
42 94.6515 84.1515
43 78.4115 94.6515
44 60.0315 78.4115
45 60.6315 60.0315
46 58.4315 60.6315
47 34.7315 58.4315
48 35.8850 34.7315
49 -17.0980 35.8850
50 -37.8780 -17.0980
51 -57.9780 -37.8780
52 -59.6485 -57.9780
53 -71.0485 -59.6485
54 -74.5485 -71.0485
55 -72.4885 -74.5485
56 -71.1685 -72.4885
57 -48.0685 -71.1685
58 -47.4685 -48.0685
59 -16.1685 -47.4685
60 1.0850 -16.1685
61 NA 1.0850
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -68.3505 -79.7675
[2,] -58.9305 -68.3505
[3,] -52.2305 -58.9305
[4,] -55.2010 -52.2305
[5,] -59.1010 -55.2010
[6,] -62.6010 -59.1010
[7,] -80.1410 -62.6010
[8,] -92.5210 -80.1410
[9,] -82.4210 -92.5210
[10,] -88.7210 -82.4210
[11,] -80.4210 -88.7210
[12,] -59.6675 -80.4210
[13,] -50.1505 -59.6675
[14,] -34.9305 -50.1505
[15,] -15.0305 -34.9305
[16,] -9.0010 -15.0305
[17,] -8.9010 -9.0010
[18,] -17.7010 -8.9010
[19,] 0.7590 -17.7010
[20,] 33.1790 0.7590
[21,] 10.9790 33.1790
[22,] 22.3790 10.9790
[23,] 27.2790 22.3790
[24,] 48.1325 27.2790
[25,] 64.5495 48.1325
[26,] 69.8695 64.5495
[27,] 77.8695 69.8695
[28,] 58.2990 77.8695
[29,] 54.8990 58.2990
[30,] 60.1990 54.8990
[31,] 73.4590 60.1990
[32,] 70.4790 73.4590
[33,] 58.8790 70.4790
[34,] 55.3790 58.8790
[35,] 34.5790 55.3790
[36,] 54.3325 34.5790
[37,] 71.0495 54.3325
[38,] 61.8695 71.0495
[39,] 47.3695 61.8695
[40,] 65.5515 47.3695
[41,] 84.1515 65.5515
[42,] 94.6515 84.1515
[43,] 78.4115 94.6515
[44,] 60.0315 78.4115
[45,] 60.6315 60.0315
[46,] 58.4315 60.6315
[47,] 34.7315 58.4315
[48,] 35.8850 34.7315
[49,] -17.0980 35.8850
[50,] -37.8780 -17.0980
[51,] -57.9780 -37.8780
[52,] -59.6485 -57.9780
[53,] -71.0485 -59.6485
[54,] -74.5485 -71.0485
[55,] -72.4885 -74.5485
[56,] -71.1685 -72.4885
[57,] -48.0685 -71.1685
[58,] -47.4685 -48.0685
[59,] -16.1685 -47.4685
[60,] 1.0850 -16.1685
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -68.3505 -79.7675
2 -58.9305 -68.3505
3 -52.2305 -58.9305
4 -55.2010 -52.2305
5 -59.1010 -55.2010
6 -62.6010 -59.1010
7 -80.1410 -62.6010
8 -92.5210 -80.1410
9 -82.4210 -92.5210
10 -88.7210 -82.4210
11 -80.4210 -88.7210
12 -59.6675 -80.4210
13 -50.1505 -59.6675
14 -34.9305 -50.1505
15 -15.0305 -34.9305
16 -9.0010 -15.0305
17 -8.9010 -9.0010
18 -17.7010 -8.9010
19 0.7590 -17.7010
20 33.1790 0.7590
21 10.9790 33.1790
22 22.3790 10.9790
23 27.2790 22.3790
24 48.1325 27.2790
25 64.5495 48.1325
26 69.8695 64.5495
27 77.8695 69.8695
28 58.2990 77.8695
29 54.8990 58.2990
30 60.1990 54.8990
31 73.4590 60.1990
32 70.4790 73.4590
33 58.8790 70.4790
34 55.3790 58.8790
35 34.5790 55.3790
36 54.3325 34.5790
37 71.0495 54.3325
38 61.8695 71.0495
39 47.3695 61.8695
40 65.5515 47.3695
41 84.1515 65.5515
42 94.6515 84.1515
43 78.4115 94.6515
44 60.0315 78.4115
45 60.6315 60.0315
46 58.4315 60.6315
47 34.7315 58.4315
48 35.8850 34.7315
49 -17.0980 35.8850
50 -37.8780 -17.0980
51 -57.9780 -37.8780
52 -59.6485 -57.9780
53 -71.0485 -59.6485
54 -74.5485 -71.0485
55 -72.4885 -74.5485
56 -71.1685 -72.4885
57 -48.0685 -71.1685
58 -47.4685 -48.0685
59 -16.1685 -47.4685
60 1.0850 -16.1685
> 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/7gw3x1259397052.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/8o9ct1259397052.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/9jumq1259397052.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/10hhhk1259397052.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/11s1341259397052.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/123qxe1259397052.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/13ubsj1259397052.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/14xlna1259397052.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/15q1581259397052.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/16fyit1259397052.tab")
+ }
>
> system("convert tmp/1rs0l1259397052.ps tmp/1rs0l1259397052.png")
> system("convert tmp/2wjbv1259397052.ps tmp/2wjbv1259397052.png")
> system("convert tmp/3f0pi1259397052.ps tmp/3f0pi1259397052.png")
> system("convert tmp/4j70k1259397052.ps tmp/4j70k1259397052.png")
> system("convert tmp/5hmon1259397052.ps tmp/5hmon1259397052.png")
> system("convert tmp/60c6i1259397052.ps tmp/60c6i1259397052.png")
> system("convert tmp/7gw3x1259397052.ps tmp/7gw3x1259397052.png")
> system("convert tmp/8o9ct1259397052.ps tmp/8o9ct1259397052.png")
> system("convert tmp/9jumq1259397052.ps tmp/9jumq1259397052.png")
> system("convert tmp/10hhhk1259397052.ps tmp/10hhhk1259397052.png")
>
>
> proc.time()
user system elapsed
2.397 1.587 7.010