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(10057,10900,11771,11992,11993,14504,11727,11477,13578,11555,11846,11397,10066,10269,14279,13870,13695,14420,11424,9704,12464,14301,13464,9893,11572,12380,16692,16052,16459,14761,13654,13480,18068,16560,14530,10650,11651,13735,13360,17818,20613,16231,13862,12004,17734,15034,12609,12320,10833,11350,13648,14890,16325,18045,15616,11926,16855,15083,12520,12355),dim=c(1,60),dimnames=list(c('Pas'),1:60))
> y <- array(NA,dim=c(1,60),dimnames=list(c('Pas'),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 = '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
Pas M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 10057 1 0 0 0 0 0 0 0 0 0 0 1
2 10900 0 1 0 0 0 0 0 0 0 0 0 2
3 11771 0 0 1 0 0 0 0 0 0 0 0 3
4 11992 0 0 0 1 0 0 0 0 0 0 0 4
5 11993 0 0 0 0 1 0 0 0 0 0 0 5
6 14504 0 0 0 0 0 1 0 0 0 0 0 6
7 11727 0 0 0 0 0 0 1 0 0 0 0 7
8 11477 0 0 0 0 0 0 0 1 0 0 0 8
9 13578 0 0 0 0 0 0 0 0 1 0 0 9
10 11555 0 0 0 0 0 0 0 0 0 1 0 10
11 11846 0 0 0 0 0 0 0 0 0 0 1 11
12 11397 0 0 0 0 0 0 0 0 0 0 0 12
13 10066 1 0 0 0 0 0 0 0 0 0 0 13
14 10269 0 1 0 0 0 0 0 0 0 0 0 14
15 14279 0 0 1 0 0 0 0 0 0 0 0 15
16 13870 0 0 0 1 0 0 0 0 0 0 0 16
17 13695 0 0 0 0 1 0 0 0 0 0 0 17
18 14420 0 0 0 0 0 1 0 0 0 0 0 18
19 11424 0 0 0 0 0 0 1 0 0 0 0 19
20 9704 0 0 0 0 0 0 0 1 0 0 0 20
21 12464 0 0 0 0 0 0 0 0 1 0 0 21
22 14301 0 0 0 0 0 0 0 0 0 1 0 22
23 13464 0 0 0 0 0 0 0 0 0 0 1 23
24 9893 0 0 0 0 0 0 0 0 0 0 0 24
25 11572 1 0 0 0 0 0 0 0 0 0 0 25
26 12380 0 1 0 0 0 0 0 0 0 0 0 26
27 16692 0 0 1 0 0 0 0 0 0 0 0 27
28 16052 0 0 0 1 0 0 0 0 0 0 0 28
29 16459 0 0 0 0 1 0 0 0 0 0 0 29
30 14761 0 0 0 0 0 1 0 0 0 0 0 30
31 13654 0 0 0 0 0 0 1 0 0 0 0 31
32 13480 0 0 0 0 0 0 0 1 0 0 0 32
33 18068 0 0 0 0 0 0 0 0 1 0 0 33
34 16560 0 0 0 0 0 0 0 0 0 1 0 34
35 14530 0 0 0 0 0 0 0 0 0 0 1 35
36 10650 0 0 0 0 0 0 0 0 0 0 0 36
37 11651 1 0 0 0 0 0 0 0 0 0 0 37
38 13735 0 1 0 0 0 0 0 0 0 0 0 38
39 13360 0 0 1 0 0 0 0 0 0 0 0 39
40 17818 0 0 0 1 0 0 0 0 0 0 0 40
41 20613 0 0 0 0 1 0 0 0 0 0 0 41
42 16231 0 0 0 0 0 1 0 0 0 0 0 42
43 13862 0 0 0 0 0 0 1 0 0 0 0 43
44 12004 0 0 0 0 0 0 0 1 0 0 0 44
45 17734 0 0 0 0 0 0 0 0 1 0 0 45
46 15034 0 0 0 0 0 0 0 0 0 1 0 46
47 12609 0 0 0 0 0 0 0 0 0 0 1 47
48 12320 0 0 0 0 0 0 0 0 0 0 0 48
49 10833 1 0 0 0 0 0 0 0 0 0 0 49
50 11350 0 1 0 0 0 0 0 0 0 0 0 50
51 13648 0 0 1 0 0 0 0 0 0 0 0 51
52 14890 0 0 0 1 0 0 0 0 0 0 0 52
53 16325 0 0 0 0 1 0 0 0 0 0 0 53
54 18045 0 0 0 0 0 1 0 0 0 0 0 54
55 15616 0 0 0 0 0 0 1 0 0 0 0 55
56 11926 0 0 0 0 0 0 0 1 0 0 0 56
57 16855 0 0 0 0 0 0 0 0 1 0 0 57
58 15083 0 0 0 0 0 0 0 0 0 1 0 58
59 12520 0 0 0 0 0 0 0 0 0 0 1 59
60 12355 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
9262.50 142.40 976.16 3142.12 4059.29 4894.65
M6 M7 M8 M9 M10 M11
4612.62 2219.78 624.14 4588.51 3298.07 1728.04
t
57.24
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2589.0 -917.2 -119.4 993.3 4109.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9262.50 785.46 11.792 1.21e-15 ***
M1 142.40 955.56 0.149 0.88218
M2 976.16 954.13 1.023 0.31150
M3 3142.12 952.83 3.298 0.00186 **
M4 4059.29 951.67 4.265 9.56e-05 ***
M5 4894.65 950.65 5.149 5.08e-06 ***
M6 4612.62 949.76 4.857 1.36e-05 ***
M7 2219.78 949.01 2.339 0.02364 *
M8 624.14 948.40 0.658 0.51368
M9 4588.51 947.92 4.841 1.44e-05 ***
M10 3298.07 947.57 3.481 0.00109 **
M11 1728.04 947.37 1.824 0.07451 .
t 57.24 11.39 5.023 7.78e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1498 on 47 degrees of freedom
Multiple R-squared: 0.6938, Adjusted R-squared: 0.6157
F-statistic: 8.876 on 12 and 47 DF, p-value: 1.587e-08
> 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.26341155 0.526823101 0.736588450
[2,] 0.16782474 0.335649483 0.832175259
[3,] 0.11039794 0.220795872 0.889602064
[4,] 0.08159705 0.163194110 0.918402945
[5,] 0.13752485 0.275049698 0.862475151
[6,] 0.23085457 0.461709146 0.769145427
[7,] 0.30059351 0.601187019 0.699406490
[8,] 0.23372310 0.467446192 0.766276904
[9,] 0.28164302 0.563286033 0.718356984
[10,] 0.21332618 0.426652370 0.786673815
[11,] 0.16272767 0.325455341 0.837272329
[12,] 0.31389287 0.627785735 0.686107132
[13,] 0.29108854 0.582177079 0.708911460
[14,] 0.35752692 0.715053834 0.642473083
[15,] 0.45994721 0.919894422 0.540052789
[16,] 0.43950739 0.879014790 0.560492605
[17,] 0.38050742 0.761014847 0.619492577
[18,] 0.46062766 0.921255314 0.539372343
[19,] 0.41128337 0.822566748 0.588716626
[20,] 0.35948949 0.718978981 0.640510510
[21,] 0.45867054 0.917341081 0.541329460
[22,] 0.38722202 0.774444039 0.612777980
[23,] 0.35415973 0.708319460 0.645840270
[24,] 0.38372417 0.767448332 0.616275834
[25,] 0.44531771 0.890635427 0.554682287
[26,] 0.95982781 0.080344378 0.040172189
[27,] 0.96848680 0.063026410 0.031513205
[28,] 0.99727807 0.005443865 0.002721933
[29,] 0.98625793 0.027484143 0.013742072
> postscript(file="/var/www/html/freestat/rcomp/tmp/1z9h01290773013.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/2aihl1290773013.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/3aihl1290773013.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/43ryo1290773013.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/53ryo1290773013.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
594.86667 546.86667 -805.33333 -1558.73333 -2450.33333 285.46667
7 8 9 10 11 12
-155.93333 1132.46667 -788.13333 -1577.93333 225.86667 1447.66667
13 14 15 16 17 18
-82.96667 -770.96667 1015.83333 -367.56667 -1435.16667 -485.36667
19 20 21 22 23 24
-1145.76667 -1327.36667 -2588.96667 481.23333 1157.03333 -743.16667
25 26 27 28 29 30
736.20000 653.20000 2742.00000 1127.60000 642.00000 -831.20000
31 32 33 34 35 36
397.40000 1761.80000 2328.20000 2053.40000 1536.20000 -673.00000
37 38 39 40 41 42
128.36667 1321.36667 -1276.83333 2206.76667 4109.16667 -48.03333
43 44 45 46 47 48
-81.43333 -401.03333 1307.36667 -159.43333 -1071.63333 310.16667
49 50 51 52 53 54
-1376.46667 -1750.46667 -1675.66667 -1408.06667 -865.66667 1079.13333
55 56 57 58 59 60
985.73333 -1165.86667 -258.46667 -797.26667 -1847.46667 -341.66667
> postscript(file="/var/www/html/freestat/rcomp/tmp/63ryo1290773013.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 594.86667 NA
1 546.86667 594.86667
2 -805.33333 546.86667
3 -1558.73333 -805.33333
4 -2450.33333 -1558.73333
5 285.46667 -2450.33333
6 -155.93333 285.46667
7 1132.46667 -155.93333
8 -788.13333 1132.46667
9 -1577.93333 -788.13333
10 225.86667 -1577.93333
11 1447.66667 225.86667
12 -82.96667 1447.66667
13 -770.96667 -82.96667
14 1015.83333 -770.96667
15 -367.56667 1015.83333
16 -1435.16667 -367.56667
17 -485.36667 -1435.16667
18 -1145.76667 -485.36667
19 -1327.36667 -1145.76667
20 -2588.96667 -1327.36667
21 481.23333 -2588.96667
22 1157.03333 481.23333
23 -743.16667 1157.03333
24 736.20000 -743.16667
25 653.20000 736.20000
26 2742.00000 653.20000
27 1127.60000 2742.00000
28 642.00000 1127.60000
29 -831.20000 642.00000
30 397.40000 -831.20000
31 1761.80000 397.40000
32 2328.20000 1761.80000
33 2053.40000 2328.20000
34 1536.20000 2053.40000
35 -673.00000 1536.20000
36 128.36667 -673.00000
37 1321.36667 128.36667
38 -1276.83333 1321.36667
39 2206.76667 -1276.83333
40 4109.16667 2206.76667
41 -48.03333 4109.16667
42 -81.43333 -48.03333
43 -401.03333 -81.43333
44 1307.36667 -401.03333
45 -159.43333 1307.36667
46 -1071.63333 -159.43333
47 310.16667 -1071.63333
48 -1376.46667 310.16667
49 -1750.46667 -1376.46667
50 -1675.66667 -1750.46667
51 -1408.06667 -1675.66667
52 -865.66667 -1408.06667
53 1079.13333 -865.66667
54 985.73333 1079.13333
55 -1165.86667 985.73333
56 -258.46667 -1165.86667
57 -797.26667 -258.46667
58 -1847.46667 -797.26667
59 -341.66667 -1847.46667
60 NA -341.66667
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 546.86667 594.86667
[2,] -805.33333 546.86667
[3,] -1558.73333 -805.33333
[4,] -2450.33333 -1558.73333
[5,] 285.46667 -2450.33333
[6,] -155.93333 285.46667
[7,] 1132.46667 -155.93333
[8,] -788.13333 1132.46667
[9,] -1577.93333 -788.13333
[10,] 225.86667 -1577.93333
[11,] 1447.66667 225.86667
[12,] -82.96667 1447.66667
[13,] -770.96667 -82.96667
[14,] 1015.83333 -770.96667
[15,] -367.56667 1015.83333
[16,] -1435.16667 -367.56667
[17,] -485.36667 -1435.16667
[18,] -1145.76667 -485.36667
[19,] -1327.36667 -1145.76667
[20,] -2588.96667 -1327.36667
[21,] 481.23333 -2588.96667
[22,] 1157.03333 481.23333
[23,] -743.16667 1157.03333
[24,] 736.20000 -743.16667
[25,] 653.20000 736.20000
[26,] 2742.00000 653.20000
[27,] 1127.60000 2742.00000
[28,] 642.00000 1127.60000
[29,] -831.20000 642.00000
[30,] 397.40000 -831.20000
[31,] 1761.80000 397.40000
[32,] 2328.20000 1761.80000
[33,] 2053.40000 2328.20000
[34,] 1536.20000 2053.40000
[35,] -673.00000 1536.20000
[36,] 128.36667 -673.00000
[37,] 1321.36667 128.36667
[38,] -1276.83333 1321.36667
[39,] 2206.76667 -1276.83333
[40,] 4109.16667 2206.76667
[41,] -48.03333 4109.16667
[42,] -81.43333 -48.03333
[43,] -401.03333 -81.43333
[44,] 1307.36667 -401.03333
[45,] -159.43333 1307.36667
[46,] -1071.63333 -159.43333
[47,] 310.16667 -1071.63333
[48,] -1376.46667 310.16667
[49,] -1750.46667 -1376.46667
[50,] -1675.66667 -1750.46667
[51,] -1408.06667 -1675.66667
[52,] -865.66667 -1408.06667
[53,] 1079.13333 -865.66667
[54,] 985.73333 1079.13333
[55,] -1165.86667 985.73333
[56,] -258.46667 -1165.86667
[57,] -797.26667 -258.46667
[58,] -1847.46667 -797.26667
[59,] -341.66667 -1847.46667
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 546.86667 594.86667
2 -805.33333 546.86667
3 -1558.73333 -805.33333
4 -2450.33333 -1558.73333
5 285.46667 -2450.33333
6 -155.93333 285.46667
7 1132.46667 -155.93333
8 -788.13333 1132.46667
9 -1577.93333 -788.13333
10 225.86667 -1577.93333
11 1447.66667 225.86667
12 -82.96667 1447.66667
13 -770.96667 -82.96667
14 1015.83333 -770.96667
15 -367.56667 1015.83333
16 -1435.16667 -367.56667
17 -485.36667 -1435.16667
18 -1145.76667 -485.36667
19 -1327.36667 -1145.76667
20 -2588.96667 -1327.36667
21 481.23333 -2588.96667
22 1157.03333 481.23333
23 -743.16667 1157.03333
24 736.20000 -743.16667
25 653.20000 736.20000
26 2742.00000 653.20000
27 1127.60000 2742.00000
28 642.00000 1127.60000
29 -831.20000 642.00000
30 397.40000 -831.20000
31 1761.80000 397.40000
32 2328.20000 1761.80000
33 2053.40000 2328.20000
34 1536.20000 2053.40000
35 -673.00000 1536.20000
36 128.36667 -673.00000
37 1321.36667 128.36667
38 -1276.83333 1321.36667
39 2206.76667 -1276.83333
40 4109.16667 2206.76667
41 -48.03333 4109.16667
42 -81.43333 -48.03333
43 -401.03333 -81.43333
44 1307.36667 -401.03333
45 -159.43333 1307.36667
46 -1071.63333 -159.43333
47 310.16667 -1071.63333
48 -1376.46667 310.16667
49 -1750.46667 -1376.46667
50 -1675.66667 -1750.46667
51 -1408.06667 -1675.66667
52 -865.66667 -1408.06667
53 1079.13333 -865.66667
54 985.73333 1079.13333
55 -1165.86667 985.73333
56 -258.46667 -1165.86667
57 -797.26667 -258.46667
58 -1847.46667 -797.26667
59 -341.66667 -1847.46667
> 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/7v1xr1290773013.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/8oaeu1290773013.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/9oaeu1290773013.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/freestat/rcomp/tmp/10oaeu1290773013.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/11savh1290773013.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/12vbb51290773013.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/132c8z1290773013.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/145c7n1290773013.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/15qdnb1290773013.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/16ud4h1290773013.tab")
+ }
>
> try(system("convert tmp/1z9h01290773013.ps tmp/1z9h01290773013.png",intern=TRUE))
character(0)
> try(system("convert tmp/2aihl1290773013.ps tmp/2aihl1290773013.png",intern=TRUE))
character(0)
> try(system("convert tmp/3aihl1290773013.ps tmp/3aihl1290773013.png",intern=TRUE))
character(0)
> try(system("convert tmp/43ryo1290773013.ps tmp/43ryo1290773013.png",intern=TRUE))
character(0)
> try(system("convert tmp/53ryo1290773013.ps tmp/53ryo1290773013.png",intern=TRUE))
character(0)
> try(system("convert tmp/63ryo1290773013.ps tmp/63ryo1290773013.png",intern=TRUE))
character(0)
> try(system("convert tmp/7v1xr1290773013.ps tmp/7v1xr1290773013.png",intern=TRUE))
character(0)
> try(system("convert tmp/8oaeu1290773013.ps tmp/8oaeu1290773013.png",intern=TRUE))
character(0)
> try(system("convert tmp/9oaeu1290773013.ps tmp/9oaeu1290773013.png",intern=TRUE))
character(0)
> try(system("convert tmp/10oaeu1290773013.ps tmp/10oaeu1290773013.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.822 2.517 4.186