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(6539,2605,6699,2682,6962,2755,6981,2760,7024,2735,6940,2659,6774,2654,6671,2670,6965,2785,6969,2845,6822,2723,6878,2746,6691,2767,6837,2940,7018,2977,7167,2993,7076,2892,7171,2824,7093,2771,6971,2686,7142,2738,7047,2723,6999,2731,6650,2632,6475,2606,6437,2605,6639,2646,6422,2627,6272,2535,6232,2456,6003,2404,5673,2319,6050,2519,5977,2504,5796,2382,5752,2394,5609,2381,5839,2501,6069,2532,6006,2515,5809,2429,5797,2389,5502,2261,5568,2272,5864,2439,5764,2373,5615,2327,5615,2364,5681,2388,5915,2553,6334,2663,6494,2694,6620,2679,6578,2611,6495,2580,6538,2627,6737,2732,6651,2707,6530,2633,6563,2683),dim=c(2,60),dimnames=list(c('Voeding-Mannen','Landbouw-Mannen'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Voeding-Mannen','Landbouw-Mannen'),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
Voeding-Mannen Landbouw-Mannen M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 6539 2605 1 0 0 0 0 0 0 0 0 0 0 1
2 6699 2682 0 1 0 0 0 0 0 0 0 0 0 2
3 6962 2755 0 0 1 0 0 0 0 0 0 0 0 3
4 6981 2760 0 0 0 1 0 0 0 0 0 0 0 4
5 7024 2735 0 0 0 0 1 0 0 0 0 0 0 5
6 6940 2659 0 0 0 0 0 1 0 0 0 0 0 6
7 6774 2654 0 0 0 0 0 0 1 0 0 0 0 7
8 6671 2670 0 0 0 0 0 0 0 1 0 0 0 8
9 6965 2785 0 0 0 0 0 0 0 0 1 0 0 9
10 6969 2845 0 0 0 0 0 0 0 0 0 1 0 10
11 6822 2723 0 0 0 0 0 0 0 0 0 0 1 11
12 6878 2746 0 0 0 0 0 0 0 0 0 0 0 12
13 6691 2767 1 0 0 0 0 0 0 0 0 0 0 13
14 6837 2940 0 1 0 0 0 0 0 0 0 0 0 14
15 7018 2977 0 0 1 0 0 0 0 0 0 0 0 15
16 7167 2993 0 0 0 1 0 0 0 0 0 0 0 16
17 7076 2892 0 0 0 0 1 0 0 0 0 0 0 17
18 7171 2824 0 0 0 0 0 1 0 0 0 0 0 18
19 7093 2771 0 0 0 0 0 0 1 0 0 0 0 19
20 6971 2686 0 0 0 0 0 0 0 1 0 0 0 20
21 7142 2738 0 0 0 0 0 0 0 0 1 0 0 21
22 7047 2723 0 0 0 0 0 0 0 0 0 1 0 22
23 6999 2731 0 0 0 0 0 0 0 0 0 0 1 23
24 6650 2632 0 0 0 0 0 0 0 0 0 0 0 24
25 6475 2606 1 0 0 0 0 0 0 0 0 0 0 25
26 6437 2605 0 1 0 0 0 0 0 0 0 0 0 26
27 6639 2646 0 0 1 0 0 0 0 0 0 0 0 27
28 6422 2627 0 0 0 1 0 0 0 0 0 0 0 28
29 6272 2535 0 0 0 0 1 0 0 0 0 0 0 29
30 6232 2456 0 0 0 0 0 1 0 0 0 0 0 30
31 6003 2404 0 0 0 0 0 0 1 0 0 0 0 31
32 5673 2319 0 0 0 0 0 0 0 1 0 0 0 32
33 6050 2519 0 0 0 0 0 0 0 0 1 0 0 33
34 5977 2504 0 0 0 0 0 0 0 0 0 1 0 34
35 5796 2382 0 0 0 0 0 0 0 0 0 0 1 35
36 5752 2394 0 0 0 0 0 0 0 0 0 0 0 36
37 5609 2381 1 0 0 0 0 0 0 0 0 0 0 37
38 5839 2501 0 1 0 0 0 0 0 0 0 0 0 38
39 6069 2532 0 0 1 0 0 0 0 0 0 0 0 39
40 6006 2515 0 0 0 1 0 0 0 0 0 0 0 40
41 5809 2429 0 0 0 0 1 0 0 0 0 0 0 41
42 5797 2389 0 0 0 0 0 1 0 0 0 0 0 42
43 5502 2261 0 0 0 0 0 0 1 0 0 0 0 43
44 5568 2272 0 0 0 0 0 0 0 1 0 0 0 44
45 5864 2439 0 0 0 0 0 0 0 0 1 0 0 45
46 5764 2373 0 0 0 0 0 0 0 0 0 1 0 46
47 5615 2327 0 0 0 0 0 0 0 0 0 0 1 47
48 5615 2364 0 0 0 0 0 0 0 0 0 0 0 48
49 5681 2388 1 0 0 0 0 0 0 0 0 0 0 49
50 5915 2553 0 1 0 0 0 0 0 0 0 0 0 50
51 6334 2663 0 0 1 0 0 0 0 0 0 0 0 51
52 6494 2694 0 0 0 1 0 0 0 0 0 0 0 52
53 6620 2679 0 0 0 0 1 0 0 0 0 0 0 53
54 6578 2611 0 0 0 0 0 1 0 0 0 0 0 54
55 6495 2580 0 0 0 0 0 0 1 0 0 0 0 55
56 6538 2627 0 0 0 0 0 0 0 1 0 0 0 56
57 6737 2732 0 0 0 0 0 0 0 0 1 0 0 57
58 6651 2707 0 0 0 0 0 0 0 0 0 1 0 58
59 6530 2633 0 0 0 0 0 0 0 0 0 0 1 59
60 6563 2683 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) `Landbouw-Mannen` M1 M2
-405.961 2.673 -101.429 -236.178
M3 M4 M5 M6
-128.966 -123.617 -2.593 162.045
M7 M8 M9 M10
139.941 106.359 36.482 3.391
M11 t
68.793 -4.302
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-339.3498 -87.9779 0.2254 90.1730 283.8007
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -405.9613 396.8517 -1.023 0.31168
`Landbouw-Mannen` 2.6728 0.1419 18.829 < 2e-16 ***
M1 -101.4290 96.8571 -1.047 0.30048
M2 -236.1780 96.4283 -2.449 0.01818 *
M3 -128.9655 97.2891 -1.326 0.19152
M4 -123.6168 97.3593 -1.270 0.21058
M5 -2.5934 96.2083 -0.027 0.97861
M6 162.0447 95.7924 1.692 0.09748 .
M7 139.9406 96.0224 1.457 0.15181
M8 106.3590 96.1448 1.106 0.27438
M9 36.4821 96.0359 0.380 0.70578
M10 3.3913 95.9085 0.035 0.97195
M11 68.7932 95.5573 0.720 0.47522
t -4.3015 1.3419 -3.206 0.00245 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 151.1 on 46 degrees of freedom
Multiple R-squared: 0.9321, Adjusted R-squared: 0.9128
F-statistic: 48.54 on 13 and 46 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.1208392 2.416784e-01 8.791608e-01
[2,] 0.1872830 3.745659e-01 8.127170e-01
[3,] 0.3146835 6.293670e-01 6.853165e-01
[4,] 0.2011976 4.023953e-01 7.988024e-01
[5,] 0.2854130 5.708260e-01 7.145870e-01
[6,] 0.4102611 8.205222e-01 5.897389e-01
[7,] 0.3065337 6.130673e-01 6.934663e-01
[8,] 0.8147368 3.705264e-01 1.852632e-01
[9,] 0.8758624 2.482752e-01 1.241376e-01
[10,] 0.9671142 6.577169e-02 3.288584e-02
[11,] 0.9934824 1.303522e-02 6.517609e-03
[12,] 0.9985097 2.980592e-03 1.490296e-03
[13,] 0.9995195 9.610054e-04 4.805027e-04
[14,] 0.9999456 1.088482e-04 5.442412e-05
[15,] 0.9999590 8.193983e-05 4.096992e-05
[16,] 0.9999579 8.412161e-05 4.206081e-05
[17,] 0.9999512 9.755975e-05 4.877988e-05
[18,] 0.9999768 4.639831e-05 2.319916e-05
[19,] 0.9999316 1.367493e-04 6.837463e-05
[20,] 0.9997985 4.030833e-04 2.015416e-04
[21,] 0.9995833 8.334909e-04 4.167455e-04
[22,] 0.9986887 2.622581e-03 1.311291e-03
[23,] 0.9991531 1.693717e-03 8.468586e-04
[24,] 0.9999170 1.659989e-04 8.299946e-05
[25,] 0.9994906 1.018813e-03 5.094066e-04
[26,] 0.9972732 5.453504e-03 2.726752e-03
[27,] 0.9956319 8.736263e-03 4.368131e-03
> postscript(file="/var/www/html/rcomp/tmp/16wrm1258725713.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/2b5y41258725713.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/3qd7y1258725713.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/49q8t1258725713.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/5fvgq1258725713.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
88.158230 181.406406 146.384145 150.973182 144.070187 102.863212
7 8 9 10 11 12
-23.367359 -131.248427 -70.437129 -189.410228 -71.434142 -3.812882
13 14 15 16 17 18
-141.210277 -318.546830 -339.349817 -234.161114 -171.934532 -55.523568
19 20 21 22 23 24
34.538226 177.605676 283.800702 266.284423 135.802021 124.499708
25 26 27 28 29 30
124.721920 228.445189 217.951171 50.686390 29.858154 40.669451
31 32 33 34 35 36
-22.941513 -87.874063 -171.247161 -166.763440 -82.787354 -85.765760
37 38 39 40 41 42
-88.289397 -39.969797 4.263761 -14.346536 -98.211317 -163.637566
43 44 45 46 47 48
-90.118953 -15.636232 -91.808330 21.986029 -65.167462 -90.964809
49 50 51 52 53 54
16.619524 -51.334968 -29.249260 46.848078 96.217508 75.628471
55 56 57 58 59 60
101.889598 57.153045 49.691919 67.903216 83.586937 56.043742
> postscript(file="/var/www/html/rcomp/tmp/6jg5f1258725713.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 88.158230 NA
1 181.406406 88.158230
2 146.384145 181.406406
3 150.973182 146.384145
4 144.070187 150.973182
5 102.863212 144.070187
6 -23.367359 102.863212
7 -131.248427 -23.367359
8 -70.437129 -131.248427
9 -189.410228 -70.437129
10 -71.434142 -189.410228
11 -3.812882 -71.434142
12 -141.210277 -3.812882
13 -318.546830 -141.210277
14 -339.349817 -318.546830
15 -234.161114 -339.349817
16 -171.934532 -234.161114
17 -55.523568 -171.934532
18 34.538226 -55.523568
19 177.605676 34.538226
20 283.800702 177.605676
21 266.284423 283.800702
22 135.802021 266.284423
23 124.499708 135.802021
24 124.721920 124.499708
25 228.445189 124.721920
26 217.951171 228.445189
27 50.686390 217.951171
28 29.858154 50.686390
29 40.669451 29.858154
30 -22.941513 40.669451
31 -87.874063 -22.941513
32 -171.247161 -87.874063
33 -166.763440 -171.247161
34 -82.787354 -166.763440
35 -85.765760 -82.787354
36 -88.289397 -85.765760
37 -39.969797 -88.289397
38 4.263761 -39.969797
39 -14.346536 4.263761
40 -98.211317 -14.346536
41 -163.637566 -98.211317
42 -90.118953 -163.637566
43 -15.636232 -90.118953
44 -91.808330 -15.636232
45 21.986029 -91.808330
46 -65.167462 21.986029
47 -90.964809 -65.167462
48 16.619524 -90.964809
49 -51.334968 16.619524
50 -29.249260 -51.334968
51 46.848078 -29.249260
52 96.217508 46.848078
53 75.628471 96.217508
54 101.889598 75.628471
55 57.153045 101.889598
56 49.691919 57.153045
57 67.903216 49.691919
58 83.586937 67.903216
59 56.043742 83.586937
60 NA 56.043742
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 181.406406 88.158230
[2,] 146.384145 181.406406
[3,] 150.973182 146.384145
[4,] 144.070187 150.973182
[5,] 102.863212 144.070187
[6,] -23.367359 102.863212
[7,] -131.248427 -23.367359
[8,] -70.437129 -131.248427
[9,] -189.410228 -70.437129
[10,] -71.434142 -189.410228
[11,] -3.812882 -71.434142
[12,] -141.210277 -3.812882
[13,] -318.546830 -141.210277
[14,] -339.349817 -318.546830
[15,] -234.161114 -339.349817
[16,] -171.934532 -234.161114
[17,] -55.523568 -171.934532
[18,] 34.538226 -55.523568
[19,] 177.605676 34.538226
[20,] 283.800702 177.605676
[21,] 266.284423 283.800702
[22,] 135.802021 266.284423
[23,] 124.499708 135.802021
[24,] 124.721920 124.499708
[25,] 228.445189 124.721920
[26,] 217.951171 228.445189
[27,] 50.686390 217.951171
[28,] 29.858154 50.686390
[29,] 40.669451 29.858154
[30,] -22.941513 40.669451
[31,] -87.874063 -22.941513
[32,] -171.247161 -87.874063
[33,] -166.763440 -171.247161
[34,] -82.787354 -166.763440
[35,] -85.765760 -82.787354
[36,] -88.289397 -85.765760
[37,] -39.969797 -88.289397
[38,] 4.263761 -39.969797
[39,] -14.346536 4.263761
[40,] -98.211317 -14.346536
[41,] -163.637566 -98.211317
[42,] -90.118953 -163.637566
[43,] -15.636232 -90.118953
[44,] -91.808330 -15.636232
[45,] 21.986029 -91.808330
[46,] -65.167462 21.986029
[47,] -90.964809 -65.167462
[48,] 16.619524 -90.964809
[49,] -51.334968 16.619524
[50,] -29.249260 -51.334968
[51,] 46.848078 -29.249260
[52,] 96.217508 46.848078
[53,] 75.628471 96.217508
[54,] 101.889598 75.628471
[55,] 57.153045 101.889598
[56,] 49.691919 57.153045
[57,] 67.903216 49.691919
[58,] 83.586937 67.903216
[59,] 56.043742 83.586937
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 181.406406 88.158230
2 146.384145 181.406406
3 150.973182 146.384145
4 144.070187 150.973182
5 102.863212 144.070187
6 -23.367359 102.863212
7 -131.248427 -23.367359
8 -70.437129 -131.248427
9 -189.410228 -70.437129
10 -71.434142 -189.410228
11 -3.812882 -71.434142
12 -141.210277 -3.812882
13 -318.546830 -141.210277
14 -339.349817 -318.546830
15 -234.161114 -339.349817
16 -171.934532 -234.161114
17 -55.523568 -171.934532
18 34.538226 -55.523568
19 177.605676 34.538226
20 283.800702 177.605676
21 266.284423 283.800702
22 135.802021 266.284423
23 124.499708 135.802021
24 124.721920 124.499708
25 228.445189 124.721920
26 217.951171 228.445189
27 50.686390 217.951171
28 29.858154 50.686390
29 40.669451 29.858154
30 -22.941513 40.669451
31 -87.874063 -22.941513
32 -171.247161 -87.874063
33 -166.763440 -171.247161
34 -82.787354 -166.763440
35 -85.765760 -82.787354
36 -88.289397 -85.765760
37 -39.969797 -88.289397
38 4.263761 -39.969797
39 -14.346536 4.263761
40 -98.211317 -14.346536
41 -163.637566 -98.211317
42 -90.118953 -163.637566
43 -15.636232 -90.118953
44 -91.808330 -15.636232
45 21.986029 -91.808330
46 -65.167462 21.986029
47 -90.964809 -65.167462
48 16.619524 -90.964809
49 -51.334968 16.619524
50 -29.249260 -51.334968
51 46.848078 -29.249260
52 96.217508 46.848078
53 75.628471 96.217508
54 101.889598 75.628471
55 57.153045 101.889598
56 49.691919 57.153045
57 67.903216 49.691919
58 83.586937 67.903216
59 56.043742 83.586937
> 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/78r411258725713.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/85a671258725713.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/9g19o1258725713.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/10ccx21258725713.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/11xz1z1258725713.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/12f7a71258725713.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/13y4rh1258725714.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/1402be1258725714.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/15j3cn1258725714.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/1643fm1258725714.tab")
+ }
>
> system("convert tmp/16wrm1258725713.ps tmp/16wrm1258725713.png")
> system("convert tmp/2b5y41258725713.ps tmp/2b5y41258725713.png")
> system("convert tmp/3qd7y1258725713.ps tmp/3qd7y1258725713.png")
> system("convert tmp/49q8t1258725713.ps tmp/49q8t1258725713.png")
> system("convert tmp/5fvgq1258725713.ps tmp/5fvgq1258725713.png")
> system("convert tmp/6jg5f1258725713.ps tmp/6jg5f1258725713.png")
> system("convert tmp/78r411258725713.ps tmp/78r411258725713.png")
> system("convert tmp/85a671258725713.ps tmp/85a671258725713.png")
> system("convert tmp/9g19o1258725713.ps tmp/9g19o1258725713.png")
> system("convert tmp/10ccx21258725713.ps tmp/10ccx21258725713.png")
>
>
> proc.time()
user system elapsed
2.411 1.576 3.842