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(9,676,8,642,9,402,9,610,9,294,9,448,10,319,9,548,9,801,9,596,8,923,9,746,9,829,9,125,9,782,9,441,9,162,9,915,10,444,10,209,9,985,9,842,9,429,10,132,9,849,9,172,10,313,9,819,9,955,10,048,10,082,10,541,10,208,10,233,9,439,9,963,10,158,9,225,10,474,9,757,10,490,10,281,10,444,10,640,10,695,10,786,9,832,9,747,10,411,9,511,10,402,9,701,10,540,10,112,10,915,11,183,10,384,10,834,9,886,10,216,10,943,9,867,10,203,10,837,10,573,10,647,11,502,10,656,10,866,10,835,9,945,10,331),dim=c(2,72),dimnames=list(c('Monthyly','births'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('Monthyly','births'),1:72))
> 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 = 'Do not include Seasonal 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
Monthyly births t
1 9 676 1
2 8 642 2
3 9 402 3
4 9 610 4
5 9 294 5
6 9 448 6
7 10 319 7
8 9 548 8
9 9 801 9
10 9 596 10
11 8 923 11
12 9 746 12
13 9 829 13
14 9 125 14
15 9 782 15
16 9 441 16
17 9 162 17
18 9 915 18
19 10 444 19
20 10 209 20
21 9 985 21
22 9 842 22
23 9 429 23
24 10 132 24
25 9 849 25
26 9 172 26
27 10 313 27
28 9 819 28
29 9 955 29
30 10 48 30
31 10 82 31
32 10 541 32
33 10 208 33
34 10 233 34
35 9 439 35
36 9 963 36
37 10 158 37
38 9 225 38
39 10 474 39
40 9 757 40
41 10 490 41
42 10 281 42
43 10 444 43
44 10 640 44
45 10 695 45
46 10 786 46
47 9 832 47
48 9 747 48
49 10 411 49
50 9 511 50
51 10 402 51
52 9 701 52
53 10 540 53
54 10 112 54
55 10 915 55
56 11 183 56
57 10 384 57
58 10 834 58
59 9 886 59
60 10 216 60
61 10 943 61
62 9 867 62
63 10 203 63
64 10 837 64
65 10 573 65
66 10 647 66
67 11 502 67
68 10 656 68
69 10 866 69
70 10 835 70
71 9 945 71
72 10 331 72
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) births t
9.403434 -0.000986 0.017940
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.86331 -0.28352 0.08328 0.24959 0.88955
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.4034341 0.1367213 68.778 < 2e-16 ***
births -0.0009860 0.0001851 -5.328 1.18e-06 ***
t 0.0179400 0.0024251 7.398 2.52e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4248 on 69 degrees of freedom
Multiple R-squared: 0.5209, Adjusted R-squared: 0.507
F-statistic: 37.5 on 2 and 69 DF, p-value: 9.467e-12
> 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.6871798 0.6256405 0.3128202
[2,] 0.7897431 0.4205139 0.2102569
[3,] 0.7170715 0.5658570 0.2829285
[4,] 0.6041982 0.7916037 0.3958018
[5,] 0.5189363 0.9621273 0.4810637
[6,] 0.6052874 0.7894251 0.3947126
[7,] 0.5148630 0.9702740 0.4851370
[8,] 0.4374739 0.8749478 0.5625261
[9,] 0.6118055 0.7763891 0.3881945
[10,] 0.5304377 0.9391247 0.4695623
[11,] 0.4614432 0.9228864 0.5385568
[12,] 0.4566600 0.9133199 0.5433400
[13,] 0.3952165 0.7904329 0.6047835
[14,] 0.5590782 0.8818435 0.4409218
[15,] 0.5492671 0.9014657 0.4507329
[16,] 0.4731767 0.9463533 0.5268233
[17,] 0.3997394 0.7994788 0.6002606
[18,] 0.4075062 0.8150125 0.5924938
[19,] 0.3604953 0.7209907 0.6395047
[20,] 0.2956661 0.5913323 0.7043339
[21,] 0.4166737 0.8333473 0.5833263
[22,] 0.4099392 0.8198783 0.5900608
[23,] 0.3494479 0.6988959 0.6505521
[24,] 0.2851485 0.5702970 0.7148515
[25,] 0.2296115 0.4592229 0.7703885
[26,] 0.1806760 0.3613521 0.8193240
[27,] 0.1963255 0.3926511 0.8036745
[28,] 0.1563206 0.3126413 0.8436794
[29,] 0.1232018 0.2464035 0.8767982
[30,] 0.1799377 0.3598753 0.8200623
[31,] 0.1425153 0.2850305 0.8574847
[32,] 0.1070486 0.2140972 0.8929514
[33,] 0.2535638 0.5071276 0.7464362
[34,] 0.2288313 0.4576626 0.7711687
[35,] 0.2253926 0.4507852 0.7746074
[36,] 0.1979867 0.3959733 0.8020133
[37,] 0.1537416 0.3074833 0.8462584
[38,] 0.1236701 0.2473402 0.8763299
[39,] 0.1164628 0.2329257 0.8835372
[40,] 0.1186950 0.2373899 0.8813050
[41,] 0.1479148 0.2958297 0.8520852
[42,] 0.1451410 0.2902820 0.8548590
[43,] 0.1591256 0.3182513 0.8408744
[44,] 0.1223568 0.2447136 0.8776432
[45,] 0.2257532 0.4515064 0.7742468
[46,] 0.1726574 0.3453149 0.8273426
[47,] 0.2636585 0.5273170 0.7363415
[48,] 0.2062575 0.4125149 0.7937425
[49,] 0.1951174 0.3902348 0.8048826
[50,] 0.1830248 0.3660497 0.8169752
[51,] 0.2747085 0.5494170 0.7252915
[52,] 0.2067956 0.4135913 0.7932044
[53,] 0.1965530 0.3931061 0.8034470
[54,] 0.2321945 0.4643891 0.7678055
[55,] 0.1875686 0.3751373 0.8124314
[56,] 0.1714962 0.3429923 0.8285038
[57,] 0.3026346 0.6052692 0.6973654
[58,] 0.4302264 0.8604528 0.5697736
[59,] 0.3163165 0.6326330 0.6836835
[60,] 0.3547991 0.7095982 0.6452009
[61,] 0.5319731 0.9360537 0.4680269
> postscript(file="/var/wessaorg/rcomp/tmp/116461322499169.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/2xu8w1322499169.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/3bign1322499169.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/4be871322499169.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/5klrw1322499169.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 = 72
Frequency = 1
1 2 3 4 5 6
0.245145625 -0.806317537 -0.060891735 0.126251270 -0.203257096 -0.069356789
7 8 9 10 11 12
0.785512338 -0.006639164 0.224872756 0.004807715 -0.690718148 0.116824136
13 14 15 16 17 18
0.180720153 -0.531346860 0.098499322 -0.255658441 -0.548685698 0.175814170
19 20 21 22 23 24
0.693479541 0.443835223 0.191012537 0.032078002 -0.393070026 0.296155151
25 26 27 28 29 30
-0.014840113 -0.700285777 0.420796842 -0.098239336 0.017913404 0.105693283
31 32 33 34 35 36
0.121276482 0.555899437 0.209629481 0.216338896 -0.598490051 -0.099778663
37 38 39 40 41 42
0.088570758 -0.863308839 0.364259177 -0.374649626 0.344154828 0.120145883
43 44 45 46 47 48
0.262919972 0.438231266 0.474519959 0.546303783 -0.426281308 -0.528029241
49 50 51 52 53 54
0.122742876 -0.796599517 0.077989128 -0.645144060 0.178173838 -0.261763828
55 56 57 58 59 60
0.512034836 0.772360497 -0.047398330 0.378350842 -0.588318393 -0.266862227
61 62 63 64 65 66
0.432002269 -0.660871881 -0.333499859 0.273668878 -0.004568741 0.050453493
67 68 69 70 71 72
0.889547006 0.023447312 0.212562269 0.164057034 -0.745425600 -0.368754783
> postscript(file="/var/wessaorg/rcomp/tmp/6ymbp1322499169.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 0.245145625 NA
1 -0.806317537 0.245145625
2 -0.060891735 -0.806317537
3 0.126251270 -0.060891735
4 -0.203257096 0.126251270
5 -0.069356789 -0.203257096
6 0.785512338 -0.069356789
7 -0.006639164 0.785512338
8 0.224872756 -0.006639164
9 0.004807715 0.224872756
10 -0.690718148 0.004807715
11 0.116824136 -0.690718148
12 0.180720153 0.116824136
13 -0.531346860 0.180720153
14 0.098499322 -0.531346860
15 -0.255658441 0.098499322
16 -0.548685698 -0.255658441
17 0.175814170 -0.548685698
18 0.693479541 0.175814170
19 0.443835223 0.693479541
20 0.191012537 0.443835223
21 0.032078002 0.191012537
22 -0.393070026 0.032078002
23 0.296155151 -0.393070026
24 -0.014840113 0.296155151
25 -0.700285777 -0.014840113
26 0.420796842 -0.700285777
27 -0.098239336 0.420796842
28 0.017913404 -0.098239336
29 0.105693283 0.017913404
30 0.121276482 0.105693283
31 0.555899437 0.121276482
32 0.209629481 0.555899437
33 0.216338896 0.209629481
34 -0.598490051 0.216338896
35 -0.099778663 -0.598490051
36 0.088570758 -0.099778663
37 -0.863308839 0.088570758
38 0.364259177 -0.863308839
39 -0.374649626 0.364259177
40 0.344154828 -0.374649626
41 0.120145883 0.344154828
42 0.262919972 0.120145883
43 0.438231266 0.262919972
44 0.474519959 0.438231266
45 0.546303783 0.474519959
46 -0.426281308 0.546303783
47 -0.528029241 -0.426281308
48 0.122742876 -0.528029241
49 -0.796599517 0.122742876
50 0.077989128 -0.796599517
51 -0.645144060 0.077989128
52 0.178173838 -0.645144060
53 -0.261763828 0.178173838
54 0.512034836 -0.261763828
55 0.772360497 0.512034836
56 -0.047398330 0.772360497
57 0.378350842 -0.047398330
58 -0.588318393 0.378350842
59 -0.266862227 -0.588318393
60 0.432002269 -0.266862227
61 -0.660871881 0.432002269
62 -0.333499859 -0.660871881
63 0.273668878 -0.333499859
64 -0.004568741 0.273668878
65 0.050453493 -0.004568741
66 0.889547006 0.050453493
67 0.023447312 0.889547006
68 0.212562269 0.023447312
69 0.164057034 0.212562269
70 -0.745425600 0.164057034
71 -0.368754783 -0.745425600
72 NA -0.368754783
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.806317537 0.245145625
[2,] -0.060891735 -0.806317537
[3,] 0.126251270 -0.060891735
[4,] -0.203257096 0.126251270
[5,] -0.069356789 -0.203257096
[6,] 0.785512338 -0.069356789
[7,] -0.006639164 0.785512338
[8,] 0.224872756 -0.006639164
[9,] 0.004807715 0.224872756
[10,] -0.690718148 0.004807715
[11,] 0.116824136 -0.690718148
[12,] 0.180720153 0.116824136
[13,] -0.531346860 0.180720153
[14,] 0.098499322 -0.531346860
[15,] -0.255658441 0.098499322
[16,] -0.548685698 -0.255658441
[17,] 0.175814170 -0.548685698
[18,] 0.693479541 0.175814170
[19,] 0.443835223 0.693479541
[20,] 0.191012537 0.443835223
[21,] 0.032078002 0.191012537
[22,] -0.393070026 0.032078002
[23,] 0.296155151 -0.393070026
[24,] -0.014840113 0.296155151
[25,] -0.700285777 -0.014840113
[26,] 0.420796842 -0.700285777
[27,] -0.098239336 0.420796842
[28,] 0.017913404 -0.098239336
[29,] 0.105693283 0.017913404
[30,] 0.121276482 0.105693283
[31,] 0.555899437 0.121276482
[32,] 0.209629481 0.555899437
[33,] 0.216338896 0.209629481
[34,] -0.598490051 0.216338896
[35,] -0.099778663 -0.598490051
[36,] 0.088570758 -0.099778663
[37,] -0.863308839 0.088570758
[38,] 0.364259177 -0.863308839
[39,] -0.374649626 0.364259177
[40,] 0.344154828 -0.374649626
[41,] 0.120145883 0.344154828
[42,] 0.262919972 0.120145883
[43,] 0.438231266 0.262919972
[44,] 0.474519959 0.438231266
[45,] 0.546303783 0.474519959
[46,] -0.426281308 0.546303783
[47,] -0.528029241 -0.426281308
[48,] 0.122742876 -0.528029241
[49,] -0.796599517 0.122742876
[50,] 0.077989128 -0.796599517
[51,] -0.645144060 0.077989128
[52,] 0.178173838 -0.645144060
[53,] -0.261763828 0.178173838
[54,] 0.512034836 -0.261763828
[55,] 0.772360497 0.512034836
[56,] -0.047398330 0.772360497
[57,] 0.378350842 -0.047398330
[58,] -0.588318393 0.378350842
[59,] -0.266862227 -0.588318393
[60,] 0.432002269 -0.266862227
[61,] -0.660871881 0.432002269
[62,] -0.333499859 -0.660871881
[63,] 0.273668878 -0.333499859
[64,] -0.004568741 0.273668878
[65,] 0.050453493 -0.004568741
[66,] 0.889547006 0.050453493
[67,] 0.023447312 0.889547006
[68,] 0.212562269 0.023447312
[69,] 0.164057034 0.212562269
[70,] -0.745425600 0.164057034
[71,] -0.368754783 -0.745425600
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.806317537 0.245145625
2 -0.060891735 -0.806317537
3 0.126251270 -0.060891735
4 -0.203257096 0.126251270
5 -0.069356789 -0.203257096
6 0.785512338 -0.069356789
7 -0.006639164 0.785512338
8 0.224872756 -0.006639164
9 0.004807715 0.224872756
10 -0.690718148 0.004807715
11 0.116824136 -0.690718148
12 0.180720153 0.116824136
13 -0.531346860 0.180720153
14 0.098499322 -0.531346860
15 -0.255658441 0.098499322
16 -0.548685698 -0.255658441
17 0.175814170 -0.548685698
18 0.693479541 0.175814170
19 0.443835223 0.693479541
20 0.191012537 0.443835223
21 0.032078002 0.191012537
22 -0.393070026 0.032078002
23 0.296155151 -0.393070026
24 -0.014840113 0.296155151
25 -0.700285777 -0.014840113
26 0.420796842 -0.700285777
27 -0.098239336 0.420796842
28 0.017913404 -0.098239336
29 0.105693283 0.017913404
30 0.121276482 0.105693283
31 0.555899437 0.121276482
32 0.209629481 0.555899437
33 0.216338896 0.209629481
34 -0.598490051 0.216338896
35 -0.099778663 -0.598490051
36 0.088570758 -0.099778663
37 -0.863308839 0.088570758
38 0.364259177 -0.863308839
39 -0.374649626 0.364259177
40 0.344154828 -0.374649626
41 0.120145883 0.344154828
42 0.262919972 0.120145883
43 0.438231266 0.262919972
44 0.474519959 0.438231266
45 0.546303783 0.474519959
46 -0.426281308 0.546303783
47 -0.528029241 -0.426281308
48 0.122742876 -0.528029241
49 -0.796599517 0.122742876
50 0.077989128 -0.796599517
51 -0.645144060 0.077989128
52 0.178173838 -0.645144060
53 -0.261763828 0.178173838
54 0.512034836 -0.261763828
55 0.772360497 0.512034836
56 -0.047398330 0.772360497
57 0.378350842 -0.047398330
58 -0.588318393 0.378350842
59 -0.266862227 -0.588318393
60 0.432002269 -0.266862227
61 -0.660871881 0.432002269
62 -0.333499859 -0.660871881
63 0.273668878 -0.333499859
64 -0.004568741 0.273668878
65 0.050453493 -0.004568741
66 0.889547006 0.050453493
67 0.023447312 0.889547006
68 0.212562269 0.023447312
69 0.164057034 0.212562269
70 -0.745425600 0.164057034
71 -0.368754783 -0.745425600
> 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/7r82n1322499169.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/8qnnr1322499169.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/9tsyr1322499169.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/10m7ox1322499169.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/11q4l11322499169.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/12jt931322499169.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/13f00d1322499169.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/14nb2e1322499169.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/15k52a1322499169.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/16sq2m1322499169.tab")
+ }
>
> try(system("convert tmp/116461322499169.ps tmp/116461322499169.png",intern=TRUE))
character(0)
> try(system("convert tmp/2xu8w1322499169.ps tmp/2xu8w1322499169.png",intern=TRUE))
character(0)
> try(system("convert tmp/3bign1322499169.ps tmp/3bign1322499169.png",intern=TRUE))
character(0)
> try(system("convert tmp/4be871322499169.ps tmp/4be871322499169.png",intern=TRUE))
character(0)
> try(system("convert tmp/5klrw1322499169.ps tmp/5klrw1322499169.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ymbp1322499169.ps tmp/6ymbp1322499169.png",intern=TRUE))
character(0)
> try(system("convert tmp/7r82n1322499169.ps tmp/7r82n1322499169.png",intern=TRUE))
character(0)
> try(system("convert tmp/8qnnr1322499169.ps tmp/8qnnr1322499169.png",intern=TRUE))
character(0)
> try(system("convert tmp/9tsyr1322499169.ps tmp/9tsyr1322499169.png",intern=TRUE))
character(0)
> try(system("convert tmp/10m7ox1322499169.ps tmp/10m7ox1322499169.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.307 0.532 3.869