R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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.
Natural language support but running in an English locale
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(62.027,56.493,65.566,62.653,53.470,59.600,42.542,42.018,44.038,44.988,43.309,26.843,69.770,64.886,79.354,63.025,54.003,55.926,45.629,40.361,43.039,44.570,43.269,25.563,68.707,60.223,74.283,61.232,61.531,65.305,51.699,44.599,35.221,55.066,45.335,28.702,69.517,69.240,71.525,77.740,62.107,65.450,51.493,43.067,49.172,54.483,38.158,27.898,58.648,56.000,62.381,59.849,48.345,55.376,45.400,38.389,44.098,48.290,41.267,31.238),dim=c(1,60),dimnames=list(c('Yt'),1:60))
> y <- array(NA,dim=c(1,60),dimnames=list(c('Yt'),1:60))
> 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
Yt M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 62.027 1 0 0 0 0 0 0 0 0 0 0
2 56.493 0 1 0 0 0 0 0 0 0 0 0
3 65.566 0 0 1 0 0 0 0 0 0 0 0
4 62.653 0 0 0 1 0 0 0 0 0 0 0
5 53.470 0 0 0 0 1 0 0 0 0 0 0
6 59.600 0 0 0 0 0 1 0 0 0 0 0
7 42.542 0 0 0 0 0 0 1 0 0 0 0
8 42.018 0 0 0 0 0 0 0 1 0 0 0
9 44.038 0 0 0 0 0 0 0 0 1 0 0
10 44.988 0 0 0 0 0 0 0 0 0 1 0
11 43.309 0 0 0 0 0 0 0 0 0 0 1
12 26.843 0 0 0 0 0 0 0 0 0 0 0
13 69.770 1 0 0 0 0 0 0 0 0 0 0
14 64.886 0 1 0 0 0 0 0 0 0 0 0
15 79.354 0 0 1 0 0 0 0 0 0 0 0
16 63.025 0 0 0 1 0 0 0 0 0 0 0
17 54.003 0 0 0 0 1 0 0 0 0 0 0
18 55.926 0 0 0 0 0 1 0 0 0 0 0
19 45.629 0 0 0 0 0 0 1 0 0 0 0
20 40.361 0 0 0 0 0 0 0 1 0 0 0
21 43.039 0 0 0 0 0 0 0 0 1 0 0
22 44.570 0 0 0 0 0 0 0 0 0 1 0
23 43.269 0 0 0 0 0 0 0 0 0 0 1
24 25.563 0 0 0 0 0 0 0 0 0 0 0
25 68.707 1 0 0 0 0 0 0 0 0 0 0
26 60.223 0 1 0 0 0 0 0 0 0 0 0
27 74.283 0 0 1 0 0 0 0 0 0 0 0
28 61.232 0 0 0 1 0 0 0 0 0 0 0
29 61.531 0 0 0 0 1 0 0 0 0 0 0
30 65.305 0 0 0 0 0 1 0 0 0 0 0
31 51.699 0 0 0 0 0 0 1 0 0 0 0
32 44.599 0 0 0 0 0 0 0 1 0 0 0
33 35.221 0 0 0 0 0 0 0 0 1 0 0
34 55.066 0 0 0 0 0 0 0 0 0 1 0
35 45.335 0 0 0 0 0 0 0 0 0 0 1
36 28.702 0 0 0 0 0 0 0 0 0 0 0
37 69.517 1 0 0 0 0 0 0 0 0 0 0
38 69.240 0 1 0 0 0 0 0 0 0 0 0
39 71.525 0 0 1 0 0 0 0 0 0 0 0
40 77.740 0 0 0 1 0 0 0 0 0 0 0
41 62.107 0 0 0 0 1 0 0 0 0 0 0
42 65.450 0 0 0 0 0 1 0 0 0 0 0
43 51.493 0 0 0 0 0 0 1 0 0 0 0
44 43.067 0 0 0 0 0 0 0 1 0 0 0
45 49.172 0 0 0 0 0 0 0 0 1 0 0
46 54.483 0 0 0 0 0 0 0 0 0 1 0
47 38.158 0 0 0 0 0 0 0 0 0 0 1
48 27.898 0 0 0 0 0 0 0 0 0 0 0
49 58.648 1 0 0 0 0 0 0 0 0 0 0
50 56.000 0 1 0 0 0 0 0 0 0 0 0
51 62.381 0 0 1 0 0 0 0 0 0 0 0
52 59.849 0 0 0 1 0 0 0 0 0 0 0
53 48.345 0 0 0 0 1 0 0 0 0 0 0
54 55.376 0 0 0 0 0 1 0 0 0 0 0
55 45.400 0 0 0 0 0 0 1 0 0 0 0
56 38.389 0 0 0 0 0 0 0 1 0 0 0
57 44.098 0 0 0 0 0 0 0 0 1 0 0
58 48.290 0 0 0 0 0 0 0 0 0 1 0
59 41.267 0 0 0 0 0 0 0 0 0 0 1
60 31.238 0 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) M1 M2 M3 M4 M5
28.05 37.68 33.32 42.57 36.85 27.84
M6 M7 M8 M9 M10 M11
32.28 19.30 13.64 15.06 21.43 14.22
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-8.2408 -3.6776 -0.4411 3.5535 12.8402
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 28.049 2.236 12.545 < 2e-16 ***
M1 37.685 3.162 11.918 5.97e-16 ***
M2 33.320 3.162 10.538 4.43e-14 ***
M3 42.573 3.162 13.464 < 2e-16 ***
M4 36.851 3.162 11.655 1.33e-15 ***
M5 27.842 3.162 8.806 1.37e-11 ***
M6 32.283 3.162 10.210 1.28e-13 ***
M7 19.304 3.162 6.105 1.73e-07 ***
M8 13.638 3.162 4.313 7.97e-05 ***
M9 15.065 3.162 4.764 1.79e-05 ***
M10 21.431 3.162 6.778 1.61e-08 ***
M11 14.219 3.162 4.497 4.37e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.999 on 48 degrees of freedom
Multiple R-squared: 0.8795, Adjusted R-squared: 0.8519
F-statistic: 31.86 on 11 and 48 DF, p-value: < 2.2e-16
> 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.84234810 0.31530380 0.1576519
[2,] 0.72899685 0.54200629 0.2710031
[3,] 0.59945958 0.80108083 0.4005404
[4,] 0.50504528 0.98990944 0.4949547
[5,] 0.40110755 0.80221510 0.5988924
[6,] 0.29262161 0.58524322 0.7073784
[7,] 0.20008855 0.40017710 0.7999115
[8,] 0.14953008 0.29906016 0.8504699
[9,] 0.09403459 0.18806917 0.9059654
[10,] 0.05907495 0.11814989 0.9409251
[11,] 0.04005316 0.08010631 0.9599468
[12,] 0.02248856 0.04497712 0.9775114
[13,] 0.01530565 0.03061129 0.9846944
[14,] 0.01047133 0.02094267 0.9895287
[15,] 0.01717449 0.03434899 0.9828255
[16,] 0.02197139 0.04394277 0.9780286
[17,] 0.02471196 0.04942392 0.9752880
[18,] 0.01667962 0.03335923 0.9833204
[19,] 0.03540542 0.07081084 0.9645946
[20,] 0.05032131 0.10064261 0.9496787
[21,] 0.03620399 0.07240799 0.9637960
[22,] 0.02163753 0.04327507 0.9783625
[23,] 0.02189720 0.04379440 0.9781028
[24,] 0.05159601 0.10319201 0.9484040
[25,] 0.04730144 0.09460289 0.9526986
[26,] 0.37623924 0.75247848 0.6237608
[27,] 0.65130085 0.69739830 0.3486992
[28,] 0.80360586 0.39278828 0.1963941
[29,] 0.81100652 0.37798696 0.1889935
[30,] 0.76449772 0.47100455 0.2355023
[31,] 0.74603383 0.50793235 0.2539662
> postscript(file="/var/www/html/freestat/rcomp/tmp/1lpqn1290957983.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/freestat/rcomp/tmp/2lpqn1290957983.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/freestat/rcomp/tmp/3vy7q1290957983.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/freestat/rcomp/tmp/4vy7q1290957983.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/freestat/rcomp/tmp/5vy7q1290957983.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 = 60
Frequency = 1
1 2 3 4 5 6 7 8 9 10
-3.7068 -4.8754 -5.0558 -2.2468 -2.4212 -0.7314 -4.8106 0.3312 0.9244 -4.4914
11 12 13 14 15 16 17 18 19 20
1.0414 -1.2058 4.0362 3.5176 8.7322 -1.8748 -1.8882 -4.4054 -1.7236 -1.3258
21 22 23 24 25 26 27 28 29 30
-0.0746 -4.9094 1.0014 -2.4858 2.9732 -1.1454 3.6612 -3.6678 5.6398 4.9736
31 32 33 34 35 36 37 38 39 40
4.3464 2.9122 -7.8926 5.5866 3.0674 0.6532 3.7832 7.8716 0.9032 12.8402
41 42 43 44 45 46 47 48 49 50
6.2158 5.1186 4.1404 1.3802 6.0584 5.0036 -4.1096 -0.1508 -7.0858 -5.3684
51 52 53 54 55 56 57 58 59 60
-8.2408 -5.0508 -7.5462 -4.9554 -1.9526 -3.2978 0.9844 -1.1894 -1.0006 3.1892
> postscript(file="/var/www/html/freestat/rcomp/tmp/6o77t1290957983.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -3.7068 NA
1 -4.8754 -3.7068
2 -5.0558 -4.8754
3 -2.2468 -5.0558
4 -2.4212 -2.2468
5 -0.7314 -2.4212
6 -4.8106 -0.7314
7 0.3312 -4.8106
8 0.9244 0.3312
9 -4.4914 0.9244
10 1.0414 -4.4914
11 -1.2058 1.0414
12 4.0362 -1.2058
13 3.5176 4.0362
14 8.7322 3.5176
15 -1.8748 8.7322
16 -1.8882 -1.8748
17 -4.4054 -1.8882
18 -1.7236 -4.4054
19 -1.3258 -1.7236
20 -0.0746 -1.3258
21 -4.9094 -0.0746
22 1.0014 -4.9094
23 -2.4858 1.0014
24 2.9732 -2.4858
25 -1.1454 2.9732
26 3.6612 -1.1454
27 -3.6678 3.6612
28 5.6398 -3.6678
29 4.9736 5.6398
30 4.3464 4.9736
31 2.9122 4.3464
32 -7.8926 2.9122
33 5.5866 -7.8926
34 3.0674 5.5866
35 0.6532 3.0674
36 3.7832 0.6532
37 7.8716 3.7832
38 0.9032 7.8716
39 12.8402 0.9032
40 6.2158 12.8402
41 5.1186 6.2158
42 4.1404 5.1186
43 1.3802 4.1404
44 6.0584 1.3802
45 5.0036 6.0584
46 -4.1096 5.0036
47 -0.1508 -4.1096
48 -7.0858 -0.1508
49 -5.3684 -7.0858
50 -8.2408 -5.3684
51 -5.0508 -8.2408
52 -7.5462 -5.0508
53 -4.9554 -7.5462
54 -1.9526 -4.9554
55 -3.2978 -1.9526
56 0.9844 -3.2978
57 -1.1894 0.9844
58 -1.0006 -1.1894
59 3.1892 -1.0006
60 NA 3.1892
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -4.8754 -3.7068
[2,] -5.0558 -4.8754
[3,] -2.2468 -5.0558
[4,] -2.4212 -2.2468
[5,] -0.7314 -2.4212
[6,] -4.8106 -0.7314
[7,] 0.3312 -4.8106
[8,] 0.9244 0.3312
[9,] -4.4914 0.9244
[10,] 1.0414 -4.4914
[11,] -1.2058 1.0414
[12,] 4.0362 -1.2058
[13,] 3.5176 4.0362
[14,] 8.7322 3.5176
[15,] -1.8748 8.7322
[16,] -1.8882 -1.8748
[17,] -4.4054 -1.8882
[18,] -1.7236 -4.4054
[19,] -1.3258 -1.7236
[20,] -0.0746 -1.3258
[21,] -4.9094 -0.0746
[22,] 1.0014 -4.9094
[23,] -2.4858 1.0014
[24,] 2.9732 -2.4858
[25,] -1.1454 2.9732
[26,] 3.6612 -1.1454
[27,] -3.6678 3.6612
[28,] 5.6398 -3.6678
[29,] 4.9736 5.6398
[30,] 4.3464 4.9736
[31,] 2.9122 4.3464
[32,] -7.8926 2.9122
[33,] 5.5866 -7.8926
[34,] 3.0674 5.5866
[35,] 0.6532 3.0674
[36,] 3.7832 0.6532
[37,] 7.8716 3.7832
[38,] 0.9032 7.8716
[39,] 12.8402 0.9032
[40,] 6.2158 12.8402
[41,] 5.1186 6.2158
[42,] 4.1404 5.1186
[43,] 1.3802 4.1404
[44,] 6.0584 1.3802
[45,] 5.0036 6.0584
[46,] -4.1096 5.0036
[47,] -0.1508 -4.1096
[48,] -7.0858 -0.1508
[49,] -5.3684 -7.0858
[50,] -8.2408 -5.3684
[51,] -5.0508 -8.2408
[52,] -7.5462 -5.0508
[53,] -4.9554 -7.5462
[54,] -1.9526 -4.9554
[55,] -3.2978 -1.9526
[56,] 0.9844 -3.2978
[57,] -1.1894 0.9844
[58,] -1.0006 -1.1894
[59,] 3.1892 -1.0006
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -4.8754 -3.7068
2 -5.0558 -4.8754
3 -2.2468 -5.0558
4 -2.4212 -2.2468
5 -0.7314 -2.4212
6 -4.8106 -0.7314
7 0.3312 -4.8106
8 0.9244 0.3312
9 -4.4914 0.9244
10 1.0414 -4.4914
11 -1.2058 1.0414
12 4.0362 -1.2058
13 3.5176 4.0362
14 8.7322 3.5176
15 -1.8748 8.7322
16 -1.8882 -1.8748
17 -4.4054 -1.8882
18 -1.7236 -4.4054
19 -1.3258 -1.7236
20 -0.0746 -1.3258
21 -4.9094 -0.0746
22 1.0014 -4.9094
23 -2.4858 1.0014
24 2.9732 -2.4858
25 -1.1454 2.9732
26 3.6612 -1.1454
27 -3.6678 3.6612
28 5.6398 -3.6678
29 4.9736 5.6398
30 4.3464 4.9736
31 2.9122 4.3464
32 -7.8926 2.9122
33 5.5866 -7.8926
34 3.0674 5.5866
35 0.6532 3.0674
36 3.7832 0.6532
37 7.8716 3.7832
38 0.9032 7.8716
39 12.8402 0.9032
40 6.2158 12.8402
41 5.1186 6.2158
42 4.1404 5.1186
43 1.3802 4.1404
44 6.0584 1.3802
45 5.0036 6.0584
46 -4.1096 5.0036
47 -0.1508 -4.1096
48 -7.0858 -0.1508
49 -5.3684 -7.0858
50 -8.2408 -5.3684
51 -5.0508 -8.2408
52 -7.5462 -5.0508
53 -4.9554 -7.5462
54 -1.9526 -4.9554
55 -3.2978 -1.9526
56 0.9844 -3.2978
57 -1.1894 0.9844
58 -1.0006 -1.1894
59 3.1892 -1.0006
> 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/freestat/rcomp/tmp/7o77t1290957983.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/freestat/rcomp/tmp/8zy6w1290957983.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/freestat/rcomp/tmp/9zy6w1290957983.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')
hat values (leverages) are all = 0.2
and there are no factor predictors; no plot no. 5
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/10s85z1290957983.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11d84n1290957983.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/freestat/rcomp/tmp/12g9kt1290957983.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/freestat/rcomp/tmp/13cii11290957983.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/freestat/rcomp/tmp/148b121290957984.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/freestat/rcomp/tmp/15uu081290957984.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/freestat/rcomp/tmp/16qlyz1290957984.tab")
+ }
> try(system("convert tmp/1lpqn1290957983.ps tmp/1lpqn1290957983.png",intern=TRUE))
character(0)
> try(system("convert tmp/2lpqn1290957983.ps tmp/2lpqn1290957983.png",intern=TRUE))
character(0)
> try(system("convert tmp/3vy7q1290957983.ps tmp/3vy7q1290957983.png",intern=TRUE))
character(0)
> try(system("convert tmp/4vy7q1290957983.ps tmp/4vy7q1290957983.png",intern=TRUE))
character(0)
> try(system("convert tmp/5vy7q1290957983.ps tmp/5vy7q1290957983.png",intern=TRUE))
character(0)
> try(system("convert tmp/6o77t1290957983.ps tmp/6o77t1290957983.png",intern=TRUE))
character(0)
> try(system("convert tmp/7o77t1290957983.ps tmp/7o77t1290957983.png",intern=TRUE))
character(0)
> try(system("convert tmp/8zy6w1290957983.ps tmp/8zy6w1290957983.png",intern=TRUE))
character(0)
> try(system("convert tmp/9zy6w1290957983.ps tmp/9zy6w1290957983.png",intern=TRUE))
character(0)
> try(system("convert tmp/10s85z1290957983.ps tmp/10s85z1290957983.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.744 2.439 4.036