R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-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(478
+ ,184
+ ,40
+ ,74
+ ,11
+ ,31
+ ,20
+ ,494
+ ,213
+ ,32
+ ,72
+ ,11
+ ,43
+ ,18
+ ,643
+ ,347
+ ,57
+ ,70
+ ,18
+ ,16
+ ,16
+ ,341
+ ,565
+ ,31
+ ,71
+ ,11
+ ,25
+ ,19
+ ,773
+ ,327
+ ,67
+ ,72
+ ,9
+ ,29
+ ,24
+ ,603
+ ,260
+ ,25
+ ,68
+ ,8
+ ,32
+ ,15
+ ,484
+ ,325
+ ,34
+ ,68
+ ,12
+ ,24
+ ,14
+ ,546
+ ,102
+ ,33
+ ,62
+ ,13
+ ,28
+ ,11
+ ,424
+ ,38
+ ,36
+ ,69
+ ,7
+ ,25
+ ,12
+ ,548
+ ,226
+ ,31
+ ,66
+ ,9
+ ,58
+ ,15
+ ,506
+ ,137
+ ,35
+ ,60
+ ,13
+ ,21
+ ,9
+ ,819
+ ,369
+ ,30
+ ,81
+ ,4
+ ,77
+ ,36
+ ,541
+ ,109
+ ,44
+ ,66
+ ,9
+ ,37
+ ,12
+ ,491
+ ,809
+ ,32
+ ,67
+ ,11
+ ,37
+ ,16
+ ,514
+ ,29
+ ,30
+ ,65
+ ,12
+ ,35
+ ,11
+ ,371
+ ,245
+ ,16
+ ,64
+ ,10
+ ,42
+ ,14
+ ,457
+ ,118
+ ,29
+ ,64
+ ,12
+ ,21
+ ,10
+ ,437
+ ,148
+ ,36
+ ,62
+ ,7
+ ,81
+ ,27
+ ,570
+ ,387
+ ,30
+ ,59
+ ,15
+ ,31
+ ,16
+ ,432
+ ,98
+ ,23
+ ,56
+ ,15
+ ,50
+ ,15
+ ,619
+ ,608
+ ,33
+ ,46
+ ,22
+ ,24
+ ,8
+ ,357
+ ,218
+ ,35
+ ,54
+ ,14
+ ,27
+ ,13
+ ,623
+ ,254
+ ,38
+ ,54
+ ,20
+ ,22
+ ,11
+ ,547
+ ,697
+ ,44
+ ,45
+ ,26
+ ,18
+ ,8
+ ,792
+ ,827
+ ,28
+ ,57
+ ,12
+ ,23
+ ,11
+ ,799
+ ,693
+ ,35
+ ,57
+ ,9
+ ,60
+ ,18
+ ,439
+ ,448
+ ,31
+ ,61
+ ,19
+ ,14
+ ,12
+ ,867
+ ,942
+ ,39
+ ,52
+ ,17
+ ,31
+ ,10
+ ,912
+ ,1017
+ ,27
+ ,44
+ ,21
+ ,24
+ ,9
+ ,462
+ ,216
+ ,36
+ ,43
+ ,18
+ ,23
+ ,8
+ ,859
+ ,673
+ ,38
+ ,48
+ ,19
+ ,22
+ ,10
+ ,805
+ ,989
+ ,46
+ ,57
+ ,14
+ ,25
+ ,12
+ ,652
+ ,630
+ ,29
+ ,47
+ ,19
+ ,25
+ ,9
+ ,776
+ ,404
+ ,32
+ ,50
+ ,19
+ ,21
+ ,9
+ ,919
+ ,692
+ ,39
+ ,48
+ ,16
+ ,32
+ ,11
+ ,732
+ ,1517
+ ,44
+ ,49
+ ,13
+ ,31
+ ,14
+ ,657
+ ,879
+ ,33
+ ,72
+ ,13
+ ,13
+ ,22
+ ,1419
+ ,631
+ ,43
+ ,59
+ ,14
+ ,21
+ ,13
+ ,989
+ ,1375
+ ,22
+ ,49
+ ,9
+ ,46
+ ,13
+ ,821
+ ,1139
+ ,30
+ ,54
+ ,13
+ ,27
+ ,12
+ ,1740
+ ,3545
+ ,86
+ ,62
+ ,22
+ ,18
+ ,15
+ ,815
+ ,706
+ ,30
+ ,47
+ ,17
+ ,39
+ ,11
+ ,760
+ ,451
+ ,32
+ ,45
+ ,34
+ ,15
+ ,10
+ ,936
+ ,433
+ ,43
+ ,48
+ ,26
+ ,23
+ ,12
+ ,863
+ ,601
+ ,20
+ ,69
+ ,23
+ ,7
+ ,12
+ ,783
+ ,1024
+ ,55
+ ,42
+ ,23
+ ,23
+ ,11
+ ,715
+ ,457
+ ,44
+ ,49
+ ,18
+ ,30
+ ,12
+ ,1504
+ ,1441
+ ,37
+ ,57
+ ,15
+ ,35
+ ,13
+ ,1324
+ ,1022
+ ,82
+ ,72
+ ,22
+ ,15
+ ,16
+ ,940
+ ,1244
+ ,66
+ ,67
+ ,26
+ ,18
+ ,16)
+ ,dim=c(7
+ ,50)
+ ,dimnames=list(c('Y1'
+ ,'X1'
+ ,'X2'
+ ,'X3'
+ ,'X4'
+ ,'X5'
+ ,'X6')
+ ,1:50))
> y <- array(NA,dim=c(7,50),dimnames=list(c('Y1','X1','X2','X3','X4','X5','X6'),1:50))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, 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
Y1 X1 X2 X3 X4 X5 X6
1 478 184 40 74 11 31 20
2 494 213 32 72 11 43 18
3 643 347 57 70 18 16 16
4 341 565 31 71 11 25 19
5 773 327 67 72 9 29 24
6 603 260 25 68 8 32 15
7 484 325 34 68 12 24 14
8 546 102 33 62 13 28 11
9 424 38 36 69 7 25 12
10 548 226 31 66 9 58 15
11 506 137 35 60 13 21 9
12 819 369 30 81 4 77 36
13 541 109 44 66 9 37 12
14 491 809 32 67 11 37 16
15 514 29 30 65 12 35 11
16 371 245 16 64 10 42 14
17 457 118 29 64 12 21 10
18 437 148 36 62 7 81 27
19 570 387 30 59 15 31 16
20 432 98 23 56 15 50 15
21 619 608 33 46 22 24 8
22 357 218 35 54 14 27 13
23 623 254 38 54 20 22 11
24 547 697 44 45 26 18 8
25 792 827 28 57 12 23 11
26 799 693 35 57 9 60 18
27 439 448 31 61 19 14 12
28 867 942 39 52 17 31 10
29 912 1017 27 44 21 24 9
30 462 216 36 43 18 23 8
31 859 673 38 48 19 22 10
32 805 989 46 57 14 25 12
33 652 630 29 47 19 25 9
34 776 404 32 50 19 21 9
35 919 692 39 48 16 32 11
36 732 1517 44 49 13 31 14
37 657 879 33 72 13 13 22
38 1419 631 43 59 14 21 13
39 989 1375 22 49 9 46 13
40 821 1139 30 54 13 27 12
41 1740 3545 86 62 22 18 15
42 815 706 30 47 17 39 11
43 760 451 32 45 34 15 10
44 936 433 43 48 26 23 12
45 863 601 20 69 23 7 12
46 783 1024 55 42 23 23 11
47 715 457 44 49 18 30 12
48 1504 1441 37 57 15 35 13
49 1324 1022 82 72 22 15 16
50 940 1244 66 67 26 18 16
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X1 X2 X3 X4 X5
100.3936 0.3323 3.9982 1.8579 7.8389 2.5588
X6
-3.2312
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-291.82 -105.31 -26.78 85.62 705.89
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 100.39361 370.69317 0.271 0.788
X1 0.33234 0.05962 5.574 1.52e-06 ***
X2 3.99817 2.68248 1.490 0.143
X3 1.85791 5.24087 0.355 0.725
X4 7.83886 7.75987 1.010 0.318
X5 2.55877 3.42695 0.747 0.459
X6 -3.23116 10.71537 -0.302 0.764
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 195.2 on 43 degrees of freedom
Multiple R-squared: 0.6132, Adjusted R-squared: 0.5592
F-statistic: 11.36 on 6 and 43 DF, p-value: 1.424e-07
> 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.173871661 0.347743323 0.8261283
[2,] 0.086401174 0.172802348 0.9135988
[3,] 0.047604626 0.095209251 0.9523954
[4,] 0.018832093 0.037664186 0.9811679
[5,] 0.009272675 0.018545349 0.9907273
[6,] 0.003411915 0.006823830 0.9965881
[7,] 0.002942144 0.005884287 0.9970579
[8,] 0.001368666 0.002737332 0.9986313
[9,] 0.010240816 0.020481632 0.9897592
[10,] 0.007752072 0.015504144 0.9922479
[11,] 0.004615962 0.009231923 0.9953840
[12,] 0.004402259 0.008804518 0.9955977
[13,] 0.004700535 0.009401070 0.9952995
[14,] 0.002829602 0.005659204 0.9971704
[15,] 0.002130114 0.004260228 0.9978699
[16,] 0.005932478 0.011864956 0.9940675
[17,] 0.004941228 0.009882456 0.9950588
[18,] 0.005944419 0.011888838 0.9940556
[19,] 0.005243569 0.010487139 0.9947564
[20,] 0.006034408 0.012068815 0.9939656
[21,] 0.004476433 0.008952865 0.9955236
[22,] 0.003798981 0.007597963 0.9962010
[23,] 0.003583901 0.007167802 0.9964161
[24,] 0.002691251 0.005382503 0.9973087
[25,] 0.003253676 0.006507352 0.9967463
[26,] 0.002143465 0.004286930 0.9978565
[27,] 0.005679242 0.011358484 0.9943208
[28,] 0.002472644 0.004945288 0.9975274
[29,] 0.276828656 0.553657311 0.7231713
[30,] 0.173496973 0.346993946 0.8265030
[31,] 0.097159517 0.194319034 0.9028405
> postscript(file="/var/wessaorg/rcomp/tmp/1kv8v1353459812.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/21kta1353459812.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/3xi451353459812.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/4h6hm1353459812.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/59x5y1353459812.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 = 50
Frequency = 1
1 2 3 4 5 6
-81.882079 -76.986177 -61.005364 -291.823520 95.078842 93.782435
7 8 9 10 11 12
-96.919451 26.569808 -41.219933 -44.558207 -17.893430 213.479251
13 14 15 16 17 18
-20.610430 -259.878862 15.178630 -134.313880 -32.951621 -146.595211
19 20 21 22 23 24
-33.776048 -94.017631 -108.874004 -192.932039 8.407685 -268.576777
25 26 27 28 29 30
81.539001 56.545123 -193.543584 20.734336 86.975268 -108.104677
31 32 33 34 35 36
120.913875 -48.830736 -44.861649 146.913213 171.761359 -275.496179
37 38 39 40 41 42
-65.310438 705.890659 106.398167 -10.415366 -167.605302 75.199880
43 44 45 46 47 48
25.583830 206.715290 195.281491 -159.240760 16.682628 505.741350
49 50
263.200669 -160.349435
> postscript(file="/var/wessaorg/rcomp/tmp/6pmhj1353459812.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 = 50
Frequency = 1
lag(myerror, k = 1) myerror
0 -81.882079 NA
1 -76.986177 -81.882079
2 -61.005364 -76.986177
3 -291.823520 -61.005364
4 95.078842 -291.823520
5 93.782435 95.078842
6 -96.919451 93.782435
7 26.569808 -96.919451
8 -41.219933 26.569808
9 -44.558207 -41.219933
10 -17.893430 -44.558207
11 213.479251 -17.893430
12 -20.610430 213.479251
13 -259.878862 -20.610430
14 15.178630 -259.878862
15 -134.313880 15.178630
16 -32.951621 -134.313880
17 -146.595211 -32.951621
18 -33.776048 -146.595211
19 -94.017631 -33.776048
20 -108.874004 -94.017631
21 -192.932039 -108.874004
22 8.407685 -192.932039
23 -268.576777 8.407685
24 81.539001 -268.576777
25 56.545123 81.539001
26 -193.543584 56.545123
27 20.734336 -193.543584
28 86.975268 20.734336
29 -108.104677 86.975268
30 120.913875 -108.104677
31 -48.830736 120.913875
32 -44.861649 -48.830736
33 146.913213 -44.861649
34 171.761359 146.913213
35 -275.496179 171.761359
36 -65.310438 -275.496179
37 705.890659 -65.310438
38 106.398167 705.890659
39 -10.415366 106.398167
40 -167.605302 -10.415366
41 75.199880 -167.605302
42 25.583830 75.199880
43 206.715290 25.583830
44 195.281491 206.715290
45 -159.240760 195.281491
46 16.682628 -159.240760
47 505.741350 16.682628
48 263.200669 505.741350
49 -160.349435 263.200669
50 NA -160.349435
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -76.986177 -81.882079
[2,] -61.005364 -76.986177
[3,] -291.823520 -61.005364
[4,] 95.078842 -291.823520
[5,] 93.782435 95.078842
[6,] -96.919451 93.782435
[7,] 26.569808 -96.919451
[8,] -41.219933 26.569808
[9,] -44.558207 -41.219933
[10,] -17.893430 -44.558207
[11,] 213.479251 -17.893430
[12,] -20.610430 213.479251
[13,] -259.878862 -20.610430
[14,] 15.178630 -259.878862
[15,] -134.313880 15.178630
[16,] -32.951621 -134.313880
[17,] -146.595211 -32.951621
[18,] -33.776048 -146.595211
[19,] -94.017631 -33.776048
[20,] -108.874004 -94.017631
[21,] -192.932039 -108.874004
[22,] 8.407685 -192.932039
[23,] -268.576777 8.407685
[24,] 81.539001 -268.576777
[25,] 56.545123 81.539001
[26,] -193.543584 56.545123
[27,] 20.734336 -193.543584
[28,] 86.975268 20.734336
[29,] -108.104677 86.975268
[30,] 120.913875 -108.104677
[31,] -48.830736 120.913875
[32,] -44.861649 -48.830736
[33,] 146.913213 -44.861649
[34,] 171.761359 146.913213
[35,] -275.496179 171.761359
[36,] -65.310438 -275.496179
[37,] 705.890659 -65.310438
[38,] 106.398167 705.890659
[39,] -10.415366 106.398167
[40,] -167.605302 -10.415366
[41,] 75.199880 -167.605302
[42,] 25.583830 75.199880
[43,] 206.715290 25.583830
[44,] 195.281491 206.715290
[45,] -159.240760 195.281491
[46,] 16.682628 -159.240760
[47,] 505.741350 16.682628
[48,] 263.200669 505.741350
[49,] -160.349435 263.200669
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -76.986177 -81.882079
2 -61.005364 -76.986177
3 -291.823520 -61.005364
4 95.078842 -291.823520
5 93.782435 95.078842
6 -96.919451 93.782435
7 26.569808 -96.919451
8 -41.219933 26.569808
9 -44.558207 -41.219933
10 -17.893430 -44.558207
11 213.479251 -17.893430
12 -20.610430 213.479251
13 -259.878862 -20.610430
14 15.178630 -259.878862
15 -134.313880 15.178630
16 -32.951621 -134.313880
17 -146.595211 -32.951621
18 -33.776048 -146.595211
19 -94.017631 -33.776048
20 -108.874004 -94.017631
21 -192.932039 -108.874004
22 8.407685 -192.932039
23 -268.576777 8.407685
24 81.539001 -268.576777
25 56.545123 81.539001
26 -193.543584 56.545123
27 20.734336 -193.543584
28 86.975268 20.734336
29 -108.104677 86.975268
30 120.913875 -108.104677
31 -48.830736 120.913875
32 -44.861649 -48.830736
33 146.913213 -44.861649
34 171.761359 146.913213
35 -275.496179 171.761359
36 -65.310438 -275.496179
37 705.890659 -65.310438
38 106.398167 705.890659
39 -10.415366 106.398167
40 -167.605302 -10.415366
41 75.199880 -167.605302
42 25.583830 75.199880
43 206.715290 25.583830
44 195.281491 206.715290
45 -159.240760 195.281491
46 16.682628 -159.240760
47 505.741350 16.682628
48 263.200669 505.741350
49 -160.349435 263.200669
> 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/7xssp1353459812.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/88aly1353459812.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/9bgi71353459812.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/10hfhg1353459812.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/11ysov1353459812.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/12lq691353459812.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/13wigw1353459812.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/14qvr51353459812.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/15j2cd1353459812.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/16fc2n1353459812.tab")
+ }
>
> try(system("convert tmp/1kv8v1353459812.ps tmp/1kv8v1353459812.png",intern=TRUE))
character(0)
> try(system("convert tmp/21kta1353459812.ps tmp/21kta1353459812.png",intern=TRUE))
character(0)
> try(system("convert tmp/3xi451353459812.ps tmp/3xi451353459812.png",intern=TRUE))
character(0)
> try(system("convert tmp/4h6hm1353459812.ps tmp/4h6hm1353459812.png",intern=TRUE))
character(0)
> try(system("convert tmp/59x5y1353459812.ps tmp/59x5y1353459812.png",intern=TRUE))
character(0)
> try(system("convert tmp/6pmhj1353459812.ps tmp/6pmhj1353459812.png",intern=TRUE))
character(0)
> try(system("convert tmp/7xssp1353459812.ps tmp/7xssp1353459812.png",intern=TRUE))
character(0)
> try(system("convert tmp/88aly1353459812.ps tmp/88aly1353459812.png",intern=TRUE))
character(0)
> try(system("convert tmp/9bgi71353459812.ps tmp/9bgi71353459812.png",intern=TRUE))
character(0)
> try(system("convert tmp/10hfhg1353459812.ps tmp/10hfhg1353459812.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
5.898 0.823 6.716