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(5560,36.68,3922,36.7,3759,36.71,4138,36.72,4634,36.73,3996,36.73,4308,36.87,4429,37.31,5219,37.39,4929,37.42,5755,37.51,5592,37.67,4163,37.67,4962,37.71,5208,37.78,4755,37.79,4491,37.84,5732,37.88,5731,38.34,5040,38.58,6102,38.72,4904,38.83,5369,38.9,5578,38.92,4619,38.94,4731,39.1,5011,39.14,5299,39.16,4146,39.32,4625,39.34,4736,39.44,4219,39.92,5116,40.19,4205,40.2,4121,40.27,5103,40.28,4300,40.3,4578,40.34,3809,40.4,5526,40.43,4247,40.48,3830,40.48,4394,40.63,4826,40.74,4409,40.77,4569,40.91,4106,40.92,4794,41.03,3914,41,3793,41.04,4405,41.33,4022,41.44,4100,41.46,4788,41.55,3163,41.55,3585,41.81,3903,41.78,4178,41.84,3863,41.84,4187,41.86),dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Y','X'),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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 5560 36.68 1 0 0 0 0 0 0 0 0 0 0 1
2 3922 36.70 0 1 0 0 0 0 0 0 0 0 0 2
3 3759 36.71 0 0 1 0 0 0 0 0 0 0 0 3
4 4138 36.72 0 0 0 1 0 0 0 0 0 0 0 4
5 4634 36.73 0 0 0 0 1 0 0 0 0 0 0 5
6 3996 36.73 0 0 0 0 0 1 0 0 0 0 0 6
7 4308 36.87 0 0 0 0 0 0 1 0 0 0 0 7
8 4429 37.31 0 0 0 0 0 0 0 1 0 0 0 8
9 5219 37.39 0 0 0 0 0 0 0 0 1 0 0 9
10 4929 37.42 0 0 0 0 0 0 0 0 0 1 0 10
11 5755 37.51 0 0 0 0 0 0 0 0 0 0 1 11
12 5592 37.67 0 0 0 0 0 0 0 0 0 0 0 12
13 4163 37.67 1 0 0 0 0 0 0 0 0 0 0 13
14 4962 37.71 0 1 0 0 0 0 0 0 0 0 0 14
15 5208 37.78 0 0 1 0 0 0 0 0 0 0 0 15
16 4755 37.79 0 0 0 1 0 0 0 0 0 0 0 16
17 4491 37.84 0 0 0 0 1 0 0 0 0 0 0 17
18 5732 37.88 0 0 0 0 0 1 0 0 0 0 0 18
19 5731 38.34 0 0 0 0 0 0 1 0 0 0 0 19
20 5040 38.58 0 0 0 0 0 0 0 1 0 0 0 20
21 6102 38.72 0 0 0 0 0 0 0 0 1 0 0 21
22 4904 38.83 0 0 0 0 0 0 0 0 0 1 0 22
23 5369 38.90 0 0 0 0 0 0 0 0 0 0 1 23
24 5578 38.92 0 0 0 0 0 0 0 0 0 0 0 24
25 4619 38.94 1 0 0 0 0 0 0 0 0 0 0 25
26 4731 39.10 0 1 0 0 0 0 0 0 0 0 0 26
27 5011 39.14 0 0 1 0 0 0 0 0 0 0 0 27
28 5299 39.16 0 0 0 1 0 0 0 0 0 0 0 28
29 4146 39.32 0 0 0 0 1 0 0 0 0 0 0 29
30 4625 39.34 0 0 0 0 0 1 0 0 0 0 0 30
31 4736 39.44 0 0 0 0 0 0 1 0 0 0 0 31
32 4219 39.92 0 0 0 0 0 0 0 1 0 0 0 32
33 5116 40.19 0 0 0 0 0 0 0 0 1 0 0 33
34 4205 40.20 0 0 0 0 0 0 0 0 0 1 0 34
35 4121 40.27 0 0 0 0 0 0 0 0 0 0 1 35
36 5103 40.28 0 0 0 0 0 0 0 0 0 0 0 36
37 4300 40.30 1 0 0 0 0 0 0 0 0 0 0 37
38 4578 40.34 0 1 0 0 0 0 0 0 0 0 0 38
39 3809 40.40 0 0 1 0 0 0 0 0 0 0 0 39
40 5526 40.43 0 0 0 1 0 0 0 0 0 0 0 40
41 4247 40.48 0 0 0 0 1 0 0 0 0 0 0 41
42 3830 40.48 0 0 0 0 0 1 0 0 0 0 0 42
43 4394 40.63 0 0 0 0 0 0 1 0 0 0 0 43
44 4826 40.74 0 0 0 0 0 0 0 1 0 0 0 44
45 4409 40.77 0 0 0 0 0 0 0 0 1 0 0 45
46 4569 40.91 0 0 0 0 0 0 0 0 0 1 0 46
47 4106 40.92 0 0 0 0 0 0 0 0 0 0 1 47
48 4794 41.03 0 0 0 0 0 0 0 0 0 0 0 48
49 3914 41.00 1 0 0 0 0 0 0 0 0 0 0 49
50 3793 41.04 0 1 0 0 0 0 0 0 0 0 0 50
51 4405 41.33 0 0 1 0 0 0 0 0 0 0 0 51
52 4022 41.44 0 0 0 1 0 0 0 0 0 0 0 52
53 4100 41.46 0 0 0 0 1 0 0 0 0 0 0 53
54 4788 41.55 0 0 0 0 0 1 0 0 0 0 0 54
55 3163 41.55 0 0 0 0 0 0 1 0 0 0 0 55
56 3585 41.81 0 0 0 0 0 0 0 1 0 0 0 56
57 3903 41.78 0 0 0 0 0 0 0 0 1 0 0 57
58 4178 41.84 0 0 0 0 0 0 0 0 0 1 0 58
59 3863 41.84 0 0 0 0 0 0 0 0 0 0 1 59
60 4187 41.86 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) X M1 M2 M3 M4
-20975.9 729.7 -739.9 -810.9 -751.5 -381.4
M5 M6 M7 M8 M9 M10
-761.3 -425.8 -590.8 -773.9 -228.6 -585.7
M11 t
-448.1 -86.8
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1039.130 -319.814 1.443 298.492 1056.439
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -20975.91 15020.83 -1.396 0.1693
X 729.66 411.03 1.775 0.0825 .
M1 -739.93 361.33 -2.048 0.0463 *
M2 -810.91 360.92 -2.247 0.0295 *
M3 -751.49 360.45 -2.085 0.0427 *
M4 -381.36 361.54 -1.055 0.2970
M5 -761.28 362.94 -2.098 0.0415 *
M6 -425.77 367.24 -1.159 0.2523
M7 -590.82 361.76 -1.633 0.1093
M8 -773.89 360.93 -2.144 0.0373 *
M9 -228.60 360.89 -0.633 0.5296
M10 -585.68 359.72 -1.628 0.1103
M11 -448.10 358.44 -1.250 0.2176
t -86.80 39.23 -2.213 0.0319 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 566.3 on 46 degrees of freedom
Multiple R-squared: 0.3819, Adjusted R-squared: 0.2073
F-statistic: 2.187 on 13 and 46 DF, p-value: 0.02597
> 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.9844791 0.03104177 0.01552088
[2,] 0.9684368 0.06312640 0.03156320
[3,] 0.9794374 0.04112516 0.02056258
[4,] 0.9636424 0.07271521 0.03635760
[5,] 0.9556008 0.08879845 0.04439922
[6,] 0.9628789 0.07424229 0.03712115
[7,] 0.9738352 0.05232962 0.02616481
[8,] 0.9596665 0.08066692 0.04033346
[9,] 0.9530616 0.09387679 0.04693840
[10,] 0.9241483 0.15170337 0.07585168
[11,] 0.8904507 0.21909863 0.10954932
[12,] 0.8402386 0.31952285 0.15976143
[13,] 0.8541302 0.29173969 0.14586984
[14,] 0.8187238 0.36255245 0.18127622
[15,] 0.7844613 0.43107733 0.21553866
[16,] 0.7671944 0.46561121 0.23280561
[17,] 0.7359221 0.52815584 0.26407792
[18,] 0.7300025 0.53999506 0.26999753
[19,] 0.7600931 0.47981375 0.23990688
[20,] 0.6758537 0.64829252 0.32414626
[21,] 0.5885785 0.82284297 0.41142148
[22,] 0.4734947 0.94698937 0.52650532
[23,] 0.5531915 0.89361710 0.44680855
[24,] 0.6468220 0.70635607 0.35317803
[25,] 0.5208725 0.95825507 0.47912754
[26,] 0.8949294 0.21014122 0.10507061
[27,] 0.9139964 0.17200727 0.08600364
> postscript(file="/var/www/html/rcomp/tmp/1p83h1258493708.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/2hx2k1258493708.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/3zpjc1258493708.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/4wj981258493708.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/5h1j81258493708.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
598.591303 -896.222172 -1039.130469 -950.759228 4.664602 -882.045504
7 8 9 10 11 12
-420.355611 -350.530469 -77.396533 54.589992 764.144141 123.096481
13 14 15 16 17 18
-479.177805 448.415457 670.727372 -72.901387 93.335918 1056.439287
19 20 21 22 23 24
971.636981 375.394747 876.748896 42.362371 405.509782 238.614960
25 26 27 28 29 30
91.747411 244.781099 522.982908 513.057518 -289.968121 -74.271490
31 32 33 34 35 36
215.604929 -381.756455 -140.258512 -614.678725 -800.531314 -187.129505
37 38 39 40 41 42
-177.997053 228.596209 -556.795245 854.982734 6.220039 -659.490067
43 44 45 46 47 48
46.903195 668.517168 -228.865740 272.857842 -248.214960 -1.779463
49 50 51 52 53 54
-33.163856 -25.570593 402.215434 -344.379637 185.747562 559.367774
55 56 57 58 59 60
-813.789494 -311.624991 -430.228111 244.868521 -120.907649 -172.802472
> postscript(file="/var/www/html/rcomp/tmp/6gj0q1258493708.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 598.591303 NA
1 -896.222172 598.591303
2 -1039.130469 -896.222172
3 -950.759228 -1039.130469
4 4.664602 -950.759228
5 -882.045504 4.664602
6 -420.355611 -882.045504
7 -350.530469 -420.355611
8 -77.396533 -350.530469
9 54.589992 -77.396533
10 764.144141 54.589992
11 123.096481 764.144141
12 -479.177805 123.096481
13 448.415457 -479.177805
14 670.727372 448.415457
15 -72.901387 670.727372
16 93.335918 -72.901387
17 1056.439287 93.335918
18 971.636981 1056.439287
19 375.394747 971.636981
20 876.748896 375.394747
21 42.362371 876.748896
22 405.509782 42.362371
23 238.614960 405.509782
24 91.747411 238.614960
25 244.781099 91.747411
26 522.982908 244.781099
27 513.057518 522.982908
28 -289.968121 513.057518
29 -74.271490 -289.968121
30 215.604929 -74.271490
31 -381.756455 215.604929
32 -140.258512 -381.756455
33 -614.678725 -140.258512
34 -800.531314 -614.678725
35 -187.129505 -800.531314
36 -177.997053 -187.129505
37 228.596209 -177.997053
38 -556.795245 228.596209
39 854.982734 -556.795245
40 6.220039 854.982734
41 -659.490067 6.220039
42 46.903195 -659.490067
43 668.517168 46.903195
44 -228.865740 668.517168
45 272.857842 -228.865740
46 -248.214960 272.857842
47 -1.779463 -248.214960
48 -33.163856 -1.779463
49 -25.570593 -33.163856
50 402.215434 -25.570593
51 -344.379637 402.215434
52 185.747562 -344.379637
53 559.367774 185.747562
54 -813.789494 559.367774
55 -311.624991 -813.789494
56 -430.228111 -311.624991
57 244.868521 -430.228111
58 -120.907649 244.868521
59 -172.802472 -120.907649
60 NA -172.802472
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -896.222172 598.591303
[2,] -1039.130469 -896.222172
[3,] -950.759228 -1039.130469
[4,] 4.664602 -950.759228
[5,] -882.045504 4.664602
[6,] -420.355611 -882.045504
[7,] -350.530469 -420.355611
[8,] -77.396533 -350.530469
[9,] 54.589992 -77.396533
[10,] 764.144141 54.589992
[11,] 123.096481 764.144141
[12,] -479.177805 123.096481
[13,] 448.415457 -479.177805
[14,] 670.727372 448.415457
[15,] -72.901387 670.727372
[16,] 93.335918 -72.901387
[17,] 1056.439287 93.335918
[18,] 971.636981 1056.439287
[19,] 375.394747 971.636981
[20,] 876.748896 375.394747
[21,] 42.362371 876.748896
[22,] 405.509782 42.362371
[23,] 238.614960 405.509782
[24,] 91.747411 238.614960
[25,] 244.781099 91.747411
[26,] 522.982908 244.781099
[27,] 513.057518 522.982908
[28,] -289.968121 513.057518
[29,] -74.271490 -289.968121
[30,] 215.604929 -74.271490
[31,] -381.756455 215.604929
[32,] -140.258512 -381.756455
[33,] -614.678725 -140.258512
[34,] -800.531314 -614.678725
[35,] -187.129505 -800.531314
[36,] -177.997053 -187.129505
[37,] 228.596209 -177.997053
[38,] -556.795245 228.596209
[39,] 854.982734 -556.795245
[40,] 6.220039 854.982734
[41,] -659.490067 6.220039
[42,] 46.903195 -659.490067
[43,] 668.517168 46.903195
[44,] -228.865740 668.517168
[45,] 272.857842 -228.865740
[46,] -248.214960 272.857842
[47,] -1.779463 -248.214960
[48,] -33.163856 -1.779463
[49,] -25.570593 -33.163856
[50,] 402.215434 -25.570593
[51,] -344.379637 402.215434
[52,] 185.747562 -344.379637
[53,] 559.367774 185.747562
[54,] -813.789494 559.367774
[55,] -311.624991 -813.789494
[56,] -430.228111 -311.624991
[57,] 244.868521 -430.228111
[58,] -120.907649 244.868521
[59,] -172.802472 -120.907649
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -896.222172 598.591303
2 -1039.130469 -896.222172
3 -950.759228 -1039.130469
4 4.664602 -950.759228
5 -882.045504 4.664602
6 -420.355611 -882.045504
7 -350.530469 -420.355611
8 -77.396533 -350.530469
9 54.589992 -77.396533
10 764.144141 54.589992
11 123.096481 764.144141
12 -479.177805 123.096481
13 448.415457 -479.177805
14 670.727372 448.415457
15 -72.901387 670.727372
16 93.335918 -72.901387
17 1056.439287 93.335918
18 971.636981 1056.439287
19 375.394747 971.636981
20 876.748896 375.394747
21 42.362371 876.748896
22 405.509782 42.362371
23 238.614960 405.509782
24 91.747411 238.614960
25 244.781099 91.747411
26 522.982908 244.781099
27 513.057518 522.982908
28 -289.968121 513.057518
29 -74.271490 -289.968121
30 215.604929 -74.271490
31 -381.756455 215.604929
32 -140.258512 -381.756455
33 -614.678725 -140.258512
34 -800.531314 -614.678725
35 -187.129505 -800.531314
36 -177.997053 -187.129505
37 228.596209 -177.997053
38 -556.795245 228.596209
39 854.982734 -556.795245
40 6.220039 854.982734
41 -659.490067 6.220039
42 46.903195 -659.490067
43 668.517168 46.903195
44 -228.865740 668.517168
45 272.857842 -228.865740
46 -248.214960 272.857842
47 -1.779463 -248.214960
48 -33.163856 -1.779463
49 -25.570593 -33.163856
50 402.215434 -25.570593
51 -344.379637 402.215434
52 185.747562 -344.379637
53 559.367774 185.747562
54 -813.789494 559.367774
55 -311.624991 -813.789494
56 -430.228111 -311.624991
57 244.868521 -430.228111
58 -120.907649 244.868521
59 -172.802472 -120.907649
> 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/7yhuv1258493708.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/8npse1258493708.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/9dllm1258493708.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/10fgdu1258493708.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/11sq6n1258493708.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/1238o51258493708.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/131ut51258493708.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/14wecj1258493709.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/15mjeh1258493709.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/16z2c91258493709.tab")
+ }
>
> system("convert tmp/1p83h1258493708.ps tmp/1p83h1258493708.png")
> system("convert tmp/2hx2k1258493708.ps tmp/2hx2k1258493708.png")
> system("convert tmp/3zpjc1258493708.ps tmp/3zpjc1258493708.png")
> system("convert tmp/4wj981258493708.ps tmp/4wj981258493708.png")
> system("convert tmp/5h1j81258493708.ps tmp/5h1j81258493708.png")
> system("convert tmp/6gj0q1258493708.ps tmp/6gj0q1258493708.png")
> system("convert tmp/7yhuv1258493708.ps tmp/7yhuv1258493708.png")
> system("convert tmp/8npse1258493708.ps tmp/8npse1258493708.png")
> system("convert tmp/9dllm1258493708.ps tmp/9dllm1258493708.png")
> system("convert tmp/10fgdu1258493708.ps tmp/10fgdu1258493708.png")
>
>
> proc.time()
user system elapsed
2.396 1.645 3.343