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(1167,333,70,669,223,44,1063,373,36,1939,873,119,678,186,30,321,111,23,2667,1277,46,345,102,39,1367,580,58,1158,420,51,1385,521,65,1155,358,40,1154,443,42,1703,690,76,1189,393,31,3083,1149,82,1357,486,36,1892,767,62,883,338,28,1627,485,38,1412,465,70,1900,816,76,777,265,33,904,307,40,2115,850,126,1858,704,56,1781,693,63,1286,387,46,1035,406,35,1557,573,108,1527,595,34,1220,394,54,1368,521,35,564,172,23,1990,835,46,1557,669,49,2057,749,56,1111,368,38,686,216,19,2011,772,29,2232,1084,26,1032,445,52,1166,451,54,1020,300,45,1735,836,56,3623,1417,596,918,330,57,1579,477,55,2790,1028,99,1496,646,51,1108,342,21,496,218,20,1750,591,58,744,255,21,1101,434,66,1612,654,47,1805,478,55,2460,753,158,1653,689,46,1234,470,45),dim=c(3,60),dimnames=list(c('TotalNrPV','TotalNrCC','TotalNrPRV'),1:60))
> y <- array(NA,dim=c(3,60),dimnames=list(c('TotalNrPV','TotalNrCC','TotalNrPRV'),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 = 'No 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
TotalNrPV TotalNrCC TotalNrPRV
1 1167 333 70
2 669 223 44
3 1063 373 36
4 1939 873 119
5 678 186 30
6 321 111 23
7 2667 1277 46
8 345 102 39
9 1367 580 58
10 1158 420 51
11 1385 521 65
12 1155 358 40
13 1154 443 42
14 1703 690 76
15 1189 393 31
16 3083 1149 82
17 1357 486 36
18 1892 767 62
19 883 338 28
20 1627 485 38
21 1412 465 70
22 1900 816 76
23 777 265 33
24 904 307 40
25 2115 850 126
26 1858 704 56
27 1781 693 63
28 1286 387 46
29 1035 406 35
30 1557 573 108
31 1527 595 34
32 1220 394 54
33 1368 521 35
34 564 172 23
35 1990 835 46
36 1557 669 49
37 2057 749 56
38 1111 368 38
39 686 216 19
40 2011 772 29
41 2232 1084 26
42 1032 445 52
43 1166 451 54
44 1020 300 45
45 1735 836 56
46 3623 1417 596
47 918 330 57
48 1579 477 55
49 2790 1028 99
50 1496 646 51
51 1108 342 21
52 496 218 20
53 1750 591 58
54 744 255 21
55 1101 434 66
56 1612 654 47
57 1805 478 55
58 2460 753 158
59 1653 689 46
60 1234 470 45
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) TotalNrCC TotalNrPRV
263.6714 2.0646 0.9807
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-309.60 -102.41 -29.94 90.64 500.51
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 263.67140 51.08374 5.162 3.23e-06 ***
TotalNrCC 2.06461 0.09828 21.008 < 2e-16 ***
TotalNrPRV 0.98068 0.36748 2.669 0.0099 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 178.2 on 57 degrees of freedom
Multiple R-squared: 0.9266, Adjusted R-squared: 0.924
F-statistic: 359.6 on 2 and 57 DF, p-value: < 2.2e-16
> 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.521375650 0.957248701 0.4786244
[2,] 0.396775338 0.793550676 0.6032247
[3,] 0.348225348 0.696450697 0.6517747
[4,] 0.236599150 0.473198300 0.7634008
[5,] 0.165009899 0.330019799 0.8349901
[6,] 0.115086507 0.230173015 0.8849135
[7,] 0.139913050 0.279826101 0.8600869
[8,] 0.087853359 0.175706718 0.9121466
[9,] 0.055624261 0.111248522 0.9443757
[10,] 0.052335802 0.104671605 0.9476642
[11,] 0.445266617 0.890533234 0.5547334
[12,] 0.375066605 0.750133210 0.6249334
[13,] 0.293633225 0.587266450 0.7063668
[14,] 0.237340919 0.474681838 0.7626591
[15,] 0.436491977 0.872983954 0.5635080
[16,] 0.395615366 0.791230732 0.6043846
[17,] 0.349149387 0.698298775 0.6508506
[18,] 0.284110670 0.568221340 0.7158893
[19,] 0.220874345 0.441748690 0.7791257
[20,] 0.166138013 0.332276026 0.8338620
[21,] 0.132040756 0.264081511 0.8679592
[22,] 0.094937706 0.189875413 0.9050623
[23,] 0.096211187 0.192422375 0.9037888
[24,] 0.075022955 0.150045910 0.9249770
[25,] 0.050977203 0.101954405 0.9490228
[26,] 0.033476248 0.066952496 0.9665238
[27,] 0.024055817 0.048111635 0.9759442
[28,] 0.014908655 0.029817310 0.9850913
[29,] 0.009864098 0.019728196 0.9901359
[30,] 0.005894597 0.011789193 0.9941054
[31,] 0.004591437 0.009182874 0.9954086
[32,] 0.004997651 0.009995303 0.9950023
[33,] 0.002901371 0.005802741 0.9970986
[34,] 0.001619733 0.003239466 0.9983803
[35,] 0.001147383 0.002294767 0.9988526
[36,] 0.003190395 0.006380791 0.9968096
[37,] 0.003642412 0.007284824 0.9963576
[38,] 0.002253881 0.004507762 0.9977461
[39,] 0.001338540 0.002677081 0.9986615
[40,] 0.009066547 0.018133093 0.9909335
[41,] 0.063269401 0.126538802 0.9367306
[42,] 0.053884810 0.107769620 0.9461152
[43,] 0.076009493 0.152018986 0.9239905
[44,] 0.079068461 0.158136921 0.9209315
[45,] 0.075663127 0.151326254 0.9243369
[46,] 0.073056855 0.146113709 0.9269431
[47,] 0.064919537 0.129839075 0.9350805
[48,] 0.050359075 0.100718150 0.9496409
[49,] 0.023448161 0.046896323 0.9765518
> postscript(file="/var/wessaorg/rcomp/tmp/1y21d1321896357.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/235wr1321896357.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/3w8hi1321896357.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/40xaa1321896357.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/5w0xd1321896357.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
147.1658969 -98.2293298 -6.0753630 -243.7767584 0.8907601 -194.3987388
7 8 9 10 11 12
-278.2894652 -167.5081368 -151.0245695 -22.8222317 -18.0773516 112.9710630
13 14 15 16 17 18
-65.4821345 -59.7838995 83.5358422 366.6761001 54.6237247 -16.0293318
19 20 21 22 23 24
-105.9685752 324.7269738 119.6373976 -122.9247398 -66.1554587 -32.7338350
25 26 27 28 29 30
-27.1554948 85.9251708 24.7711162 178.2132951 -101.2268074 4.3936790
31 32 33 34 35 36
1.5426127 89.9155830 -5.6569394 -77.3399393 -42.7319145 -135.9487218
37 38 39 40 41 42
192.0177279 50.2863254 -42.2600507 125.0100725 -295.2061573 -201.4181582
43 44 45 46 47 48
-81.7671781 92.8150319 -309.6033285 -150.7090675 -82.8914283 276.5722856
49 50 51 52 53 54
306.8223241 -151.4240563 117.6377482 -237.3699508 209.2647222 -66.7411954
55 56 57 58 59 60
-123.4369757 -48.0182134 500.5076757 486.7298869 -78.2988775 -44.1686414
> postscript(file="/var/wessaorg/rcomp/tmp/66nqk1321896357.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 147.1658969 NA
1 -98.2293298 147.1658969
2 -6.0753630 -98.2293298
3 -243.7767584 -6.0753630
4 0.8907601 -243.7767584
5 -194.3987388 0.8907601
6 -278.2894652 -194.3987388
7 -167.5081368 -278.2894652
8 -151.0245695 -167.5081368
9 -22.8222317 -151.0245695
10 -18.0773516 -22.8222317
11 112.9710630 -18.0773516
12 -65.4821345 112.9710630
13 -59.7838995 -65.4821345
14 83.5358422 -59.7838995
15 366.6761001 83.5358422
16 54.6237247 366.6761001
17 -16.0293318 54.6237247
18 -105.9685752 -16.0293318
19 324.7269738 -105.9685752
20 119.6373976 324.7269738
21 -122.9247398 119.6373976
22 -66.1554587 -122.9247398
23 -32.7338350 -66.1554587
24 -27.1554948 -32.7338350
25 85.9251708 -27.1554948
26 24.7711162 85.9251708
27 178.2132951 24.7711162
28 -101.2268074 178.2132951
29 4.3936790 -101.2268074
30 1.5426127 4.3936790
31 89.9155830 1.5426127
32 -5.6569394 89.9155830
33 -77.3399393 -5.6569394
34 -42.7319145 -77.3399393
35 -135.9487218 -42.7319145
36 192.0177279 -135.9487218
37 50.2863254 192.0177279
38 -42.2600507 50.2863254
39 125.0100725 -42.2600507
40 -295.2061573 125.0100725
41 -201.4181582 -295.2061573
42 -81.7671781 -201.4181582
43 92.8150319 -81.7671781
44 -309.6033285 92.8150319
45 -150.7090675 -309.6033285
46 -82.8914283 -150.7090675
47 276.5722856 -82.8914283
48 306.8223241 276.5722856
49 -151.4240563 306.8223241
50 117.6377482 -151.4240563
51 -237.3699508 117.6377482
52 209.2647222 -237.3699508
53 -66.7411954 209.2647222
54 -123.4369757 -66.7411954
55 -48.0182134 -123.4369757
56 500.5076757 -48.0182134
57 486.7298869 500.5076757
58 -78.2988775 486.7298869
59 -44.1686414 -78.2988775
60 NA -44.1686414
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -98.2293298 147.1658969
[2,] -6.0753630 -98.2293298
[3,] -243.7767584 -6.0753630
[4,] 0.8907601 -243.7767584
[5,] -194.3987388 0.8907601
[6,] -278.2894652 -194.3987388
[7,] -167.5081368 -278.2894652
[8,] -151.0245695 -167.5081368
[9,] -22.8222317 -151.0245695
[10,] -18.0773516 -22.8222317
[11,] 112.9710630 -18.0773516
[12,] -65.4821345 112.9710630
[13,] -59.7838995 -65.4821345
[14,] 83.5358422 -59.7838995
[15,] 366.6761001 83.5358422
[16,] 54.6237247 366.6761001
[17,] -16.0293318 54.6237247
[18,] -105.9685752 -16.0293318
[19,] 324.7269738 -105.9685752
[20,] 119.6373976 324.7269738
[21,] -122.9247398 119.6373976
[22,] -66.1554587 -122.9247398
[23,] -32.7338350 -66.1554587
[24,] -27.1554948 -32.7338350
[25,] 85.9251708 -27.1554948
[26,] 24.7711162 85.9251708
[27,] 178.2132951 24.7711162
[28,] -101.2268074 178.2132951
[29,] 4.3936790 -101.2268074
[30,] 1.5426127 4.3936790
[31,] 89.9155830 1.5426127
[32,] -5.6569394 89.9155830
[33,] -77.3399393 -5.6569394
[34,] -42.7319145 -77.3399393
[35,] -135.9487218 -42.7319145
[36,] 192.0177279 -135.9487218
[37,] 50.2863254 192.0177279
[38,] -42.2600507 50.2863254
[39,] 125.0100725 -42.2600507
[40,] -295.2061573 125.0100725
[41,] -201.4181582 -295.2061573
[42,] -81.7671781 -201.4181582
[43,] 92.8150319 -81.7671781
[44,] -309.6033285 92.8150319
[45,] -150.7090675 -309.6033285
[46,] -82.8914283 -150.7090675
[47,] 276.5722856 -82.8914283
[48,] 306.8223241 276.5722856
[49,] -151.4240563 306.8223241
[50,] 117.6377482 -151.4240563
[51,] -237.3699508 117.6377482
[52,] 209.2647222 -237.3699508
[53,] -66.7411954 209.2647222
[54,] -123.4369757 -66.7411954
[55,] -48.0182134 -123.4369757
[56,] 500.5076757 -48.0182134
[57,] 486.7298869 500.5076757
[58,] -78.2988775 486.7298869
[59,] -44.1686414 -78.2988775
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -98.2293298 147.1658969
2 -6.0753630 -98.2293298
3 -243.7767584 -6.0753630
4 0.8907601 -243.7767584
5 -194.3987388 0.8907601
6 -278.2894652 -194.3987388
7 -167.5081368 -278.2894652
8 -151.0245695 -167.5081368
9 -22.8222317 -151.0245695
10 -18.0773516 -22.8222317
11 112.9710630 -18.0773516
12 -65.4821345 112.9710630
13 -59.7838995 -65.4821345
14 83.5358422 -59.7838995
15 366.6761001 83.5358422
16 54.6237247 366.6761001
17 -16.0293318 54.6237247
18 -105.9685752 -16.0293318
19 324.7269738 -105.9685752
20 119.6373976 324.7269738
21 -122.9247398 119.6373976
22 -66.1554587 -122.9247398
23 -32.7338350 -66.1554587
24 -27.1554948 -32.7338350
25 85.9251708 -27.1554948
26 24.7711162 85.9251708
27 178.2132951 24.7711162
28 -101.2268074 178.2132951
29 4.3936790 -101.2268074
30 1.5426127 4.3936790
31 89.9155830 1.5426127
32 -5.6569394 89.9155830
33 -77.3399393 -5.6569394
34 -42.7319145 -77.3399393
35 -135.9487218 -42.7319145
36 192.0177279 -135.9487218
37 50.2863254 192.0177279
38 -42.2600507 50.2863254
39 125.0100725 -42.2600507
40 -295.2061573 125.0100725
41 -201.4181582 -295.2061573
42 -81.7671781 -201.4181582
43 92.8150319 -81.7671781
44 -309.6033285 92.8150319
45 -150.7090675 -309.6033285
46 -82.8914283 -150.7090675
47 276.5722856 -82.8914283
48 306.8223241 276.5722856
49 -151.4240563 306.8223241
50 117.6377482 -151.4240563
51 -237.3699508 117.6377482
52 209.2647222 -237.3699508
53 -66.7411954 209.2647222
54 -123.4369757 -66.7411954
55 -48.0182134 -123.4369757
56 500.5076757 -48.0182134
57 486.7298869 500.5076757
58 -78.2988775 486.7298869
59 -44.1686414 -78.2988775
> 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/7gskf1321896357.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/8rpeq1321896357.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/9sg3w1321896357.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/10mzoi1321896357.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/1186191321896357.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/123bjz1321896357.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/13bhdi1321896357.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/14442l1321896357.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/15wn0o1321896357.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/16cfe51321896357.tab")
+ }
>
> try(system("convert tmp/1y21d1321896357.ps tmp/1y21d1321896357.png",intern=TRUE))
character(0)
> try(system("convert tmp/235wr1321896357.ps tmp/235wr1321896357.png",intern=TRUE))
character(0)
> try(system("convert tmp/3w8hi1321896357.ps tmp/3w8hi1321896357.png",intern=TRUE))
character(0)
> try(system("convert tmp/40xaa1321896357.ps tmp/40xaa1321896357.png",intern=TRUE))
character(0)
> try(system("convert tmp/5w0xd1321896357.ps tmp/5w0xd1321896357.png",intern=TRUE))
character(0)
> try(system("convert tmp/66nqk1321896357.ps tmp/66nqk1321896357.png",intern=TRUE))
character(0)
> try(system("convert tmp/7gskf1321896357.ps tmp/7gskf1321896357.png",intern=TRUE))
character(0)
> try(system("convert tmp/8rpeq1321896357.ps tmp/8rpeq1321896357.png",intern=TRUE))
character(0)
> try(system("convert tmp/9sg3w1321896357.ps tmp/9sg3w1321896357.png",intern=TRUE))
character(0)
> try(system("convert tmp/10mzoi1321896357.ps tmp/10mzoi1321896357.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.326 0.534 3.925