R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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(135094,135411,135698,135880,135891,135971,136173,136358,136514,136506,136711,136891,137094,137182,137400,137479,137620,137687,137638,137612,137681,137772,137899,137983,137996,137913,137841,137656,137423,137245,137014,136747,136313,135804,135002,134383,133563,132837,132041,131381,130995,130493,130193,129962,129726,129505,129450,129320,129281,129246,129438,129715,130173,129981,129932,129873,129844,130015,130108,130260),dim=c(1,60),dimnames=list(c('y'),1:60))
> y <- array(NA,dim=c(1,60),dimnames=list(c('y'),1:60))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par20 = ''
> par19 = ''
> par18 = ''
> par17 = ''
> par16 = ''
> par15 = ''
> par14 = ''
> par13 = ''
> par12 = ''
> par11 = ''
> par10 = ''
> par9 = ''
> par8 = ''
> par7 = ''
> par6 = ''
> par5 = ''
> par4 = ''
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> library(lattice)
> library(lmtest)
Loading required package: zoo
> 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 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 135094 1 0 0 0 0 0 0 0 0 0 0 1
2 135411 0 1 0 0 0 0 0 0 0 0 0 2
3 135698 0 0 1 0 0 0 0 0 0 0 0 3
4 135880 0 0 0 1 0 0 0 0 0 0 0 4
5 135891 0 0 0 0 1 0 0 0 0 0 0 5
6 135971 0 0 0 0 0 1 0 0 0 0 0 6
7 136173 0 0 0 0 0 0 1 0 0 0 0 7
8 136358 0 0 0 0 0 0 0 1 0 0 0 8
9 136514 0 0 0 0 0 0 0 0 1 0 0 9
10 136506 0 0 0 0 0 0 0 0 0 1 0 10
11 136711 0 0 0 0 0 0 0 0 0 0 1 11
12 136891 0 0 0 0 0 0 0 0 0 0 0 12
13 137094 1 0 0 0 0 0 0 0 0 0 0 13
14 137182 0 1 0 0 0 0 0 0 0 0 0 14
15 137400 0 0 1 0 0 0 0 0 0 0 0 15
16 137479 0 0 0 1 0 0 0 0 0 0 0 16
17 137620 0 0 0 0 1 0 0 0 0 0 0 17
18 137687 0 0 0 0 0 1 0 0 0 0 0 18
19 137638 0 0 0 0 0 0 1 0 0 0 0 19
20 137612 0 0 0 0 0 0 0 1 0 0 0 20
21 137681 0 0 0 0 0 0 0 0 1 0 0 21
22 137772 0 0 0 0 0 0 0 0 0 1 0 22
23 137899 0 0 0 0 0 0 0 0 0 0 1 23
24 137983 0 0 0 0 0 0 0 0 0 0 0 24
25 137996 1 0 0 0 0 0 0 0 0 0 0 25
26 137913 0 1 0 0 0 0 0 0 0 0 0 26
27 137841 0 0 1 0 0 0 0 0 0 0 0 27
28 137656 0 0 0 1 0 0 0 0 0 0 0 28
29 137423 0 0 0 0 1 0 0 0 0 0 0 29
30 137245 0 0 0 0 0 1 0 0 0 0 0 30
31 137014 0 0 0 0 0 0 1 0 0 0 0 31
32 136747 0 0 0 0 0 0 0 1 0 0 0 32
33 136313 0 0 0 0 0 0 0 0 1 0 0 33
34 135804 0 0 0 0 0 0 0 0 0 1 0 34
35 135002 0 0 0 0 0 0 0 0 0 0 1 35
36 134383 0 0 0 0 0 0 0 0 0 0 0 36
37 133563 1 0 0 0 0 0 0 0 0 0 0 37
38 132837 0 1 0 0 0 0 0 0 0 0 0 38
39 132041 0 0 1 0 0 0 0 0 0 0 0 39
40 131381 0 0 0 1 0 0 0 0 0 0 0 40
41 130995 0 0 0 0 1 0 0 0 0 0 0 41
42 130493 0 0 0 0 0 1 0 0 0 0 0 42
43 130193 0 0 0 0 0 0 1 0 0 0 0 43
44 129962 0 0 0 0 0 0 0 1 0 0 0 44
45 129726 0 0 0 0 0 0 0 0 1 0 0 45
46 129505 0 0 0 0 0 0 0 0 0 1 0 46
47 129450 0 0 0 0 0 0 0 0 0 0 1 47
48 129320 0 0 0 0 0 0 0 0 0 0 0 48
49 129281 1 0 0 0 0 0 0 0 0 0 0 49
50 129246 0 1 0 0 0 0 0 0 0 0 0 50
51 129438 0 0 1 0 0 0 0 0 0 0 0 51
52 129715 0 0 0 1 0 0 0 0 0 0 0 52
53 130173 0 0 0 0 1 0 0 0 0 0 0 53
54 129981 0 0 0 0 0 1 0 0 0 0 0 54
55 129932 0 0 0 0 0 0 1 0 0 0 0 55
56 129873 0 0 0 0 0 0 0 1 0 0 0 56
57 129844 0 0 0 0 0 0 0 0 1 0 0 57
58 130015 0 0 0 0 0 0 0 0 0 1 0 58
59 130108 0 0 0 0 0 0 0 0 0 0 1 59
60 130260 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
139568.53 -934.37 -861.02 -734.08 -634.34 -475.00
M6 M7 M8 M9 M10 M11
-458.85 -383.11 -301.57 -235.23 -169.28 -94.54
t
-161.14
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3379.0 -1523.7 -337.1 1527.7 3395.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 139568.53 1132.34 123.256 < 2e-16 ***
M1 -934.37 1377.56 -0.678 0.501
M2 -861.02 1375.50 -0.626 0.534
M3 -734.08 1373.64 -0.534 0.596
M4 -634.34 1371.96 -0.462 0.646
M5 -475.00 1370.49 -0.347 0.730
M6 -458.85 1369.21 -0.335 0.739
M7 -383.11 1368.12 -0.280 0.781
M8 -301.57 1367.24 -0.221 0.826
M9 -235.23 1366.55 -0.172 0.864
M10 -169.28 1366.05 -0.124 0.902
M11 -94.54 1365.76 -0.069 0.945
t -161.14 16.43 -9.810 5.93e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2159 on 47 degrees of freedom
Multiple R-squared: 0.674, Adjusted R-squared: 0.5908
F-statistic: 8.099 on 12 and 47 DF, p-value: 6.036e-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,] 6.074673e-04 1.214935e-03 0.9993925327
[2,] 3.937124e-05 7.874248e-05 0.9999606288
[3,] 2.345772e-06 4.691544e-06 0.9999976542
[4,] 7.079703e-07 1.415941e-06 0.9999992920
[5,] 6.533320e-07 1.306664e-06 0.9999993467
[6,] 3.816204e-07 7.632408e-07 0.9999996184
[7,] 8.218494e-08 1.643699e-07 0.9999999178
[8,] 2.062192e-08 4.124384e-08 0.9999999794
[9,] 6.476950e-09 1.295390e-08 0.9999999935
[10,] 1.592665e-09 3.185330e-09 0.9999999984
[11,] 1.230928e-09 2.461855e-09 0.9999999988
[12,] 3.163122e-09 6.326245e-09 0.9999999968
[13,] 1.480808e-08 2.961616e-08 0.9999999852
[14,] 8.069530e-08 1.613906e-07 0.9999999193
[15,] 4.497429e-07 8.994858e-07 0.9999995503
[16,] 3.169257e-06 6.338514e-06 0.9999968307
[17,] 3.011336e-05 6.022672e-05 0.9999698866
[18,] 4.745294e-04 9.490588e-04 0.9995254706
[19,] 7.121477e-03 1.424295e-02 0.9928785229
[20,] 7.887752e-02 1.577550e-01 0.9211224827
[21,] 3.596043e-01 7.192086e-01 0.6403956917
[22,] 7.770673e-01 4.458654e-01 0.2229327064
[23,] 9.699357e-01 6.012856e-02 0.0300642806
[24,] 9.973882e-01 5.223596e-03 0.0026117982
[25,] 9.996166e-01 7.667735e-04 0.0003833868
[26,] 9.997172e-01 5.656680e-04 0.0002828340
[27,] 9.996402e-01 7.195507e-04 0.0003597753
[28,] 9.992617e-01 1.476659e-03 0.0007383293
[29,] 9.981832e-01 3.633522e-03 0.0018167611
> postscript(file="/var/wessaorg/rcomp/tmp/1hi3z1322596130.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/2rmz61322596130.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/3a2am1322596130.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/4niar1322596130.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/5t7o51322596130.ps",horizontal=F,onefile=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
-3379.01667 -2974.21667 -2653.01667 -2409.61667 -2396.81667 -2171.81667
7 8 9 10 11 12
-1884.41667 -1619.81667 -1369.01667 -1281.81667 -990.41667 -743.81667
13 14 15 16 17 18
554.69167 730.49167 982.69167 1123.09167 1265.89167 1477.89167
19 20 21 22 23 24
1514.29167 1567.89167 1731.69167 1917.89167 2131.29167 2281.89167
25 26 27 28 29 30
3390.40000 3395.20000 3357.40000 3233.80000 3002.60000 2969.60000
31 32 33 34 35 36
2824.00000 2636.60000 2297.40000 1883.60000 1168.00000 615.60000
37 38 39 40 41 42
891.10833 252.90833 -508.89167 -1107.49167 -1491.69167 -1848.69167
43 44 45 46 47 48
-2063.29167 -2214.69167 -2355.89167 -2481.69167 -2450.29167 -2513.69167
49 50 51 52 53 54
-1457.18333 -1404.38333 -1178.18333 -839.78333 -379.98333 -426.98333
55 56 57 58 59 60
-390.58333 -369.98333 -304.18333 -37.98333 141.41667 360.01667
> postscript(file="/var/wessaorg/rcomp/tmp/6fzaf1322596130.ps",horizontal=F,onefile=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 -3379.01667 NA
1 -2974.21667 -3379.01667
2 -2653.01667 -2974.21667
3 -2409.61667 -2653.01667
4 -2396.81667 -2409.61667
5 -2171.81667 -2396.81667
6 -1884.41667 -2171.81667
7 -1619.81667 -1884.41667
8 -1369.01667 -1619.81667
9 -1281.81667 -1369.01667
10 -990.41667 -1281.81667
11 -743.81667 -990.41667
12 554.69167 -743.81667
13 730.49167 554.69167
14 982.69167 730.49167
15 1123.09167 982.69167
16 1265.89167 1123.09167
17 1477.89167 1265.89167
18 1514.29167 1477.89167
19 1567.89167 1514.29167
20 1731.69167 1567.89167
21 1917.89167 1731.69167
22 2131.29167 1917.89167
23 2281.89167 2131.29167
24 3390.40000 2281.89167
25 3395.20000 3390.40000
26 3357.40000 3395.20000
27 3233.80000 3357.40000
28 3002.60000 3233.80000
29 2969.60000 3002.60000
30 2824.00000 2969.60000
31 2636.60000 2824.00000
32 2297.40000 2636.60000
33 1883.60000 2297.40000
34 1168.00000 1883.60000
35 615.60000 1168.00000
36 891.10833 615.60000
37 252.90833 891.10833
38 -508.89167 252.90833
39 -1107.49167 -508.89167
40 -1491.69167 -1107.49167
41 -1848.69167 -1491.69167
42 -2063.29167 -1848.69167
43 -2214.69167 -2063.29167
44 -2355.89167 -2214.69167
45 -2481.69167 -2355.89167
46 -2450.29167 -2481.69167
47 -2513.69167 -2450.29167
48 -1457.18333 -2513.69167
49 -1404.38333 -1457.18333
50 -1178.18333 -1404.38333
51 -839.78333 -1178.18333
52 -379.98333 -839.78333
53 -426.98333 -379.98333
54 -390.58333 -426.98333
55 -369.98333 -390.58333
56 -304.18333 -369.98333
57 -37.98333 -304.18333
58 141.41667 -37.98333
59 360.01667 141.41667
60 NA 360.01667
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2974.21667 -3379.01667
[2,] -2653.01667 -2974.21667
[3,] -2409.61667 -2653.01667
[4,] -2396.81667 -2409.61667
[5,] -2171.81667 -2396.81667
[6,] -1884.41667 -2171.81667
[7,] -1619.81667 -1884.41667
[8,] -1369.01667 -1619.81667
[9,] -1281.81667 -1369.01667
[10,] -990.41667 -1281.81667
[11,] -743.81667 -990.41667
[12,] 554.69167 -743.81667
[13,] 730.49167 554.69167
[14,] 982.69167 730.49167
[15,] 1123.09167 982.69167
[16,] 1265.89167 1123.09167
[17,] 1477.89167 1265.89167
[18,] 1514.29167 1477.89167
[19,] 1567.89167 1514.29167
[20,] 1731.69167 1567.89167
[21,] 1917.89167 1731.69167
[22,] 2131.29167 1917.89167
[23,] 2281.89167 2131.29167
[24,] 3390.40000 2281.89167
[25,] 3395.20000 3390.40000
[26,] 3357.40000 3395.20000
[27,] 3233.80000 3357.40000
[28,] 3002.60000 3233.80000
[29,] 2969.60000 3002.60000
[30,] 2824.00000 2969.60000
[31,] 2636.60000 2824.00000
[32,] 2297.40000 2636.60000
[33,] 1883.60000 2297.40000
[34,] 1168.00000 1883.60000
[35,] 615.60000 1168.00000
[36,] 891.10833 615.60000
[37,] 252.90833 891.10833
[38,] -508.89167 252.90833
[39,] -1107.49167 -508.89167
[40,] -1491.69167 -1107.49167
[41,] -1848.69167 -1491.69167
[42,] -2063.29167 -1848.69167
[43,] -2214.69167 -2063.29167
[44,] -2355.89167 -2214.69167
[45,] -2481.69167 -2355.89167
[46,] -2450.29167 -2481.69167
[47,] -2513.69167 -2450.29167
[48,] -1457.18333 -2513.69167
[49,] -1404.38333 -1457.18333
[50,] -1178.18333 -1404.38333
[51,] -839.78333 -1178.18333
[52,] -379.98333 -839.78333
[53,] -426.98333 -379.98333
[54,] -390.58333 -426.98333
[55,] -369.98333 -390.58333
[56,] -304.18333 -369.98333
[57,] -37.98333 -304.18333
[58,] 141.41667 -37.98333
[59,] 360.01667 141.41667
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2974.21667 -3379.01667
2 -2653.01667 -2974.21667
3 -2409.61667 -2653.01667
4 -2396.81667 -2409.61667
5 -2171.81667 -2396.81667
6 -1884.41667 -2171.81667
7 -1619.81667 -1884.41667
8 -1369.01667 -1619.81667
9 -1281.81667 -1369.01667
10 -990.41667 -1281.81667
11 -743.81667 -990.41667
12 554.69167 -743.81667
13 730.49167 554.69167
14 982.69167 730.49167
15 1123.09167 982.69167
16 1265.89167 1123.09167
17 1477.89167 1265.89167
18 1514.29167 1477.89167
19 1567.89167 1514.29167
20 1731.69167 1567.89167
21 1917.89167 1731.69167
22 2131.29167 1917.89167
23 2281.89167 2131.29167
24 3390.40000 2281.89167
25 3395.20000 3390.40000
26 3357.40000 3395.20000
27 3233.80000 3357.40000
28 3002.60000 3233.80000
29 2969.60000 3002.60000
30 2824.00000 2969.60000
31 2636.60000 2824.00000
32 2297.40000 2636.60000
33 1883.60000 2297.40000
34 1168.00000 1883.60000
35 615.60000 1168.00000
36 891.10833 615.60000
37 252.90833 891.10833
38 -508.89167 252.90833
39 -1107.49167 -508.89167
40 -1491.69167 -1107.49167
41 -1848.69167 -1491.69167
42 -2063.29167 -1848.69167
43 -2214.69167 -2063.29167
44 -2355.89167 -2214.69167
45 -2481.69167 -2355.89167
46 -2450.29167 -2481.69167
47 -2513.69167 -2450.29167
48 -1457.18333 -2513.69167
49 -1404.38333 -1457.18333
50 -1178.18333 -1404.38333
51 -839.78333 -1178.18333
52 -379.98333 -839.78333
53 -426.98333 -379.98333
54 -390.58333 -426.98333
55 -369.98333 -390.58333
56 -304.18333 -369.98333
57 -37.98333 -304.18333
58 141.41667 -37.98333
59 360.01667 141.41667
> 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/wessaorg/rcomp/tmp/7cse61322596130.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/88q3p1322596130.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/986p01322596130.ps",horizontal=F,onefile=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/wessaorg/rcomp/tmp/10to631322596130.ps",horizontal=F,onefile=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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11qhpe1322596130.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/wessaorg/rcomp/tmp/12005u1322596130.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/wessaorg/rcomp/tmp/137tzn1322596130.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/wessaorg/rcomp/tmp/14ta351322596130.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/wessaorg/rcomp/tmp/152hcb1322596130.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/wessaorg/rcomp/tmp/167f0g1322596130.tab")
+ }
>
> try(system("convert tmp/1hi3z1322596130.ps tmp/1hi3z1322596130.png",intern=TRUE))
character(0)
> try(system("convert tmp/2rmz61322596130.ps tmp/2rmz61322596130.png",intern=TRUE))
character(0)
> try(system("convert tmp/3a2am1322596130.ps tmp/3a2am1322596130.png",intern=TRUE))
character(0)
> try(system("convert tmp/4niar1322596130.ps tmp/4niar1322596130.png",intern=TRUE))
character(0)
> try(system("convert tmp/5t7o51322596130.ps tmp/5t7o51322596130.png",intern=TRUE))
character(0)
> try(system("convert tmp/6fzaf1322596130.ps tmp/6fzaf1322596130.png",intern=TRUE))
character(0)
> try(system("convert tmp/7cse61322596130.ps tmp/7cse61322596130.png",intern=TRUE))
character(0)
> try(system("convert tmp/88q3p1322596130.ps tmp/88q3p1322596130.png",intern=TRUE))
character(0)
> try(system("convert tmp/986p01322596130.ps tmp/986p01322596130.png",intern=TRUE))
character(0)
> try(system("convert tmp/10to631322596130.ps tmp/10to631322596130.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.162 0.537 3.821