R version 2.12.0 (2010-10-15)
Copyright (C) 2010 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(105.31
+ ,1576.23
+ ,29.29
+ ,105.63
+ ,1546.37
+ ,28.99
+ ,106.02
+ ,1545.05
+ ,28.91
+ ,105.85
+ ,1552.34
+ ,29.29
+ ,106.57
+ ,1594.3
+ ,30.96
+ ,106.48
+ ,1605.78
+ ,30.57
+ ,106.60
+ ,1673.21
+ ,30.59
+ ,106.75
+ ,1612.94
+ ,31.39
+ ,106.69
+ ,1566.34
+ ,31.28
+ ,106.69
+ ,1530.17
+ ,31.1
+ ,106.93
+ ,1582.54
+ ,31.7
+ ,107.21
+ ,1702.16
+ ,32.57
+ ,107.88
+ ,1701.93
+ ,32.49
+ ,108.84
+ ,1811.15
+ ,32.46
+ ,108.96
+ ,1924.2
+ ,32.3
+ ,109.52
+ ,2034.25
+ ,32.97
+ ,108.45
+ ,2011.13
+ ,32.9
+ ,108.67
+ ,2013.04
+ ,32.93
+ ,108.96
+ ,2151.67
+ ,33.72
+ ,108.76
+ ,1902.09
+ ,33.33
+ ,107.85
+ ,1944.01
+ ,33.44
+ ,108.78
+ ,1916.67
+ ,33.89
+ ,107.51
+ ,1967.31
+ ,34.34
+ ,108.83
+ ,2119.88
+ ,33.56
+ ,111.54
+ ,2216.38
+ ,32.67
+ ,111.74
+ ,2522.83
+ ,32.57
+ ,112.04
+ ,2647.64
+ ,33.23
+ ,111.74
+ ,2631.23
+ ,32.85
+ ,111.81
+ ,2693.41
+ ,32.61
+ ,111.86
+ ,3021.76
+ ,32.57
+ ,114.23
+ ,2953.67
+ ,32.98
+ ,114.80
+ ,2796.8
+ ,31.33
+ ,115.17
+ ,2672.05
+ ,29.8
+ ,115.11
+ ,2251.23
+ ,28.06
+ ,114.43
+ ,2046.08
+ ,25.47
+ ,114.66
+ ,2420.04
+ ,24.65
+ ,115.11
+ ,2608.89
+ ,23.94
+ ,117.74
+ ,2660.47
+ ,23.89
+ ,118.18
+ ,2493.98
+ ,23.54
+ ,118.56
+ ,2541.7
+ ,24.28
+ ,117.63
+ ,2554.6
+ ,25.51
+ ,117.71
+ ,2699.61
+ ,27.03
+ ,117.46
+ ,2805.48
+ ,27.09
+ ,117.37
+ ,2956.66
+ ,27.3
+ ,117.34
+ ,3149.51
+ ,27.11
+ ,117.09
+ ,3372.5
+ ,26.39
+ ,116.65
+ ,3379.33
+ ,27.54
+ ,116.71
+ ,3517.54
+ ,26.85
+ ,116.82
+ ,3527.34
+ ,26.82
+ ,117.33
+ ,3281.06
+ ,25.9
+ ,117.95
+ ,3089.65
+ ,24.96
+ ,123.53
+ ,3222.76
+ ,25.4
+ ,124.91
+ ,3165.76
+ ,24.38
+ ,125.99
+ ,3232.43
+ ,24.73
+ ,126.29
+ ,3229.54
+ ,25.43
+ ,125.68
+ ,3071.74
+ ,26.04
+ ,125.52
+ ,2850.17
+ ,25.59)
+ ,dim=c(3
+ ,57)
+ ,dimnames=list(c('PC&S'
+ ,'PCacao'
+ ,'PSuiker')
+ ,1:57))
> y <- array(NA,dim=c(3,57),dimnames=list(c('PC&S','PCacao','PSuiker'),1:57))
> 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'
> #'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
> 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
PC&S PCacao PSuiker t
1 105.31 1576.23 29.29 1
2 105.63 1546.37 28.99 2
3 106.02 1545.05 28.91 3
4 105.85 1552.34 29.29 4
5 106.57 1594.30 30.96 5
6 106.48 1605.78 30.57 6
7 106.60 1673.21 30.59 7
8 106.75 1612.94 31.39 8
9 106.69 1566.34 31.28 9
10 106.69 1530.17 31.10 10
11 106.93 1582.54 31.70 11
12 107.21 1702.16 32.57 12
13 107.88 1701.93 32.49 13
14 108.84 1811.15 32.46 14
15 108.96 1924.20 32.30 15
16 109.52 2034.25 32.97 16
17 108.45 2011.13 32.90 17
18 108.67 2013.04 32.93 18
19 108.96 2151.67 33.72 19
20 108.76 1902.09 33.33 20
21 107.85 1944.01 33.44 21
22 108.78 1916.67 33.89 22
23 107.51 1967.31 34.34 23
24 108.83 2119.88 33.56 24
25 111.54 2216.38 32.67 25
26 111.74 2522.83 32.57 26
27 112.04 2647.64 33.23 27
28 111.74 2631.23 32.85 28
29 111.81 2693.41 32.61 29
30 111.86 3021.76 32.57 30
31 114.23 2953.67 32.98 31
32 114.80 2796.80 31.33 32
33 115.17 2672.05 29.80 33
34 115.11 2251.23 28.06 34
35 114.43 2046.08 25.47 35
36 114.66 2420.04 24.65 36
37 115.11 2608.89 23.94 37
38 117.74 2660.47 23.89 38
39 118.18 2493.98 23.54 39
40 118.56 2541.70 24.28 40
41 117.63 2554.60 25.51 41
42 117.71 2699.61 27.03 42
43 117.46 2805.48 27.09 43
44 117.37 2956.66 27.30 44
45 117.34 3149.51 27.11 45
46 117.09 3372.50 26.39 46
47 116.65 3379.33 27.54 47
48 116.71 3517.54 26.85 48
49 116.82 3527.34 26.82 49
50 117.33 3281.06 25.90 50
51 117.95 3089.65 24.96 51
52 123.53 3222.76 25.40 52
53 124.91 3165.76 24.38 53
54 125.99 3232.43 24.73 54
55 126.29 3229.54 25.43 55
56 125.68 3071.74 26.04 56
57 125.52 2850.17 25.59 57
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) PCacao PSuiker t
114.487751 -0.001273 -0.283151 0.346452
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-3.2070 -0.8685 0.0469 0.7754 4.0583
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.145e+02 3.094e+00 37.007 < 2e-16 ***
PCacao -1.273e-03 9.985e-04 -1.275 0.20799
PSuiker -2.832e-01 9.561e-02 -2.962 0.00457 **
t 3.465e-01 4.333e-02 7.995 1.14e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.685 on 53 degrees of freedom
Multiple R-squared: 0.9245, Adjusted R-squared: 0.9203
F-statistic: 216.5 on 3 and 53 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,] 1.903666e-03 3.807333e-03 0.9980963
[2,] 3.966527e-04 7.933054e-04 0.9996033
[3,] 9.132870e-05 1.826574e-04 0.9999087
[4,] 1.244275e-05 2.488550e-05 0.9999876
[5,] 1.516897e-06 3.033795e-06 0.9999985
[6,] 2.327447e-07 4.654893e-07 0.9999998
[7,] 1.099956e-07 2.199913e-07 0.9999999
[8,] 1.524792e-07 3.049585e-07 0.9999998
[9,] 2.929209e-08 5.858419e-08 1.0000000
[10,] 5.485974e-09 1.097195e-08 1.0000000
[11,] 1.308764e-07 2.617529e-07 0.9999999
[12,] 9.157912e-08 1.831582e-07 0.9999999
[13,] 8.284141e-08 1.656828e-07 0.9999999
[14,] 2.043516e-08 4.087032e-08 1.0000000
[15,] 1.080681e-07 2.161362e-07 0.9999999
[16,] 2.887353e-08 5.774706e-08 1.0000000
[17,] 9.252100e-07 1.850420e-06 0.9999991
[18,] 6.122737e-07 1.224547e-06 0.9999994
[19,] 7.063094e-06 1.412619e-05 0.9999929
[20,] 2.475650e-06 4.951301e-06 0.9999975
[21,] 8.583463e-07 1.716693e-06 0.9999991
[22,] 3.192367e-07 6.384733e-07 0.9999997
[23,] 1.325770e-07 2.651541e-07 0.9999999
[24,] 2.458243e-07 4.916486e-07 0.9999998
[25,] 6.256171e-07 1.251234e-06 0.9999994
[26,] 4.738065e-06 9.476130e-06 0.9999953
[27,] 3.545998e-05 7.091996e-05 0.9999645
[28,] 3.560538e-05 7.121076e-05 0.9999644
[29,] 4.511648e-05 9.023296e-05 0.9999549
[30,] 7.884885e-05 1.576977e-04 0.9999212
[31,] 1.026758e-04 2.053516e-04 0.9998973
[32,] 8.857727e-05 1.771545e-04 0.9999114
[33,] 7.557293e-05 1.511459e-04 0.9999244
[34,] 6.851952e-05 1.370390e-04 0.9999315
[35,] 3.105641e-05 6.211283e-05 0.9999689
[36,] 1.618719e-05 3.237438e-05 0.9999838
[37,] 7.685664e-06 1.537133e-05 0.9999923
[38,] 8.403101e-06 1.680620e-05 0.9999916
[39,] 1.547922e-04 3.095843e-04 0.9998452
[40,] 2.245990e-03 4.491980e-03 0.9977540
[41,] 2.096274e-01 4.192549e-01 0.7903726
[42,] 2.942569e-01 5.885137e-01 0.7057431
[43,] 2.383806e-01 4.767613e-01 0.7616194
[44,] 2.060258e-01 4.120515e-01 0.7939742
> postscript(file="/var/www/rcomp/tmp/1cpon1292931240.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/www/rcomp/tmp/2cpon1292931240.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/www/rcomp/tmp/34gn81292931240.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/www/rcomp/tmp/44gn81292931240.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/www/rcomp/tmp/54gn81292931240.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 = 57
Frequency = 1
1 2 3 4 5 6
0.77542748 0.62602641 0.64524269 0.24566673 1.14548149 0.61321216
7 8 9 10 11 12
0.47824484 0.43160550 -0.06530269 -0.50875669 -0.37866407 -0.04652853
13 14 15 16 17 18
0.25407505 0.99813825 0.87026645 1.41359163 -0.05210644 -0.16763258
19 20 21 22 23 24
0.17604579 -0.79848666 -1.97043816 -1.29426877 -2.71885063 -1.77197703
25 26 27 28 29 30
0.46238709 0.67765345 0.97693286 0.20199815 -0.06327032 0.04685824
31 32 33 34 35 36
2.09983719 1.65653095 1.08808342 -0.34664775 -2.36756400 -2.24024316
37 38 39 40 41 42
-2.09737377 0.23766532 0.02021138 0.32402686 -0.58773077 -0.23923214
43 44 45 46 47 48
-0.68394907 -0.86852523 -1.05332641 -1.56983698 -2.02197225 -2.32789171
49 50 51 52 53 54
-2.56036494 -2.97076728 -3.20699706 2.32055295 2.99274091 3.91024609
55 56 57
4.05832184 3.07375291 2.15788103
> postscript(file="/var/www/rcomp/tmp/6f84t1292931240.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 = 57
Frequency = 1
lag(myerror, k = 1) myerror
0 0.77542748 NA
1 0.62602641 0.77542748
2 0.64524269 0.62602641
3 0.24566673 0.64524269
4 1.14548149 0.24566673
5 0.61321216 1.14548149
6 0.47824484 0.61321216
7 0.43160550 0.47824484
8 -0.06530269 0.43160550
9 -0.50875669 -0.06530269
10 -0.37866407 -0.50875669
11 -0.04652853 -0.37866407
12 0.25407505 -0.04652853
13 0.99813825 0.25407505
14 0.87026645 0.99813825
15 1.41359163 0.87026645
16 -0.05210644 1.41359163
17 -0.16763258 -0.05210644
18 0.17604579 -0.16763258
19 -0.79848666 0.17604579
20 -1.97043816 -0.79848666
21 -1.29426877 -1.97043816
22 -2.71885063 -1.29426877
23 -1.77197703 -2.71885063
24 0.46238709 -1.77197703
25 0.67765345 0.46238709
26 0.97693286 0.67765345
27 0.20199815 0.97693286
28 -0.06327032 0.20199815
29 0.04685824 -0.06327032
30 2.09983719 0.04685824
31 1.65653095 2.09983719
32 1.08808342 1.65653095
33 -0.34664775 1.08808342
34 -2.36756400 -0.34664775
35 -2.24024316 -2.36756400
36 -2.09737377 -2.24024316
37 0.23766532 -2.09737377
38 0.02021138 0.23766532
39 0.32402686 0.02021138
40 -0.58773077 0.32402686
41 -0.23923214 -0.58773077
42 -0.68394907 -0.23923214
43 -0.86852523 -0.68394907
44 -1.05332641 -0.86852523
45 -1.56983698 -1.05332641
46 -2.02197225 -1.56983698
47 -2.32789171 -2.02197225
48 -2.56036494 -2.32789171
49 -2.97076728 -2.56036494
50 -3.20699706 -2.97076728
51 2.32055295 -3.20699706
52 2.99274091 2.32055295
53 3.91024609 2.99274091
54 4.05832184 3.91024609
55 3.07375291 4.05832184
56 2.15788103 3.07375291
57 NA 2.15788103
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.62602641 0.77542748
[2,] 0.64524269 0.62602641
[3,] 0.24566673 0.64524269
[4,] 1.14548149 0.24566673
[5,] 0.61321216 1.14548149
[6,] 0.47824484 0.61321216
[7,] 0.43160550 0.47824484
[8,] -0.06530269 0.43160550
[9,] -0.50875669 -0.06530269
[10,] -0.37866407 -0.50875669
[11,] -0.04652853 -0.37866407
[12,] 0.25407505 -0.04652853
[13,] 0.99813825 0.25407505
[14,] 0.87026645 0.99813825
[15,] 1.41359163 0.87026645
[16,] -0.05210644 1.41359163
[17,] -0.16763258 -0.05210644
[18,] 0.17604579 -0.16763258
[19,] -0.79848666 0.17604579
[20,] -1.97043816 -0.79848666
[21,] -1.29426877 -1.97043816
[22,] -2.71885063 -1.29426877
[23,] -1.77197703 -2.71885063
[24,] 0.46238709 -1.77197703
[25,] 0.67765345 0.46238709
[26,] 0.97693286 0.67765345
[27,] 0.20199815 0.97693286
[28,] -0.06327032 0.20199815
[29,] 0.04685824 -0.06327032
[30,] 2.09983719 0.04685824
[31,] 1.65653095 2.09983719
[32,] 1.08808342 1.65653095
[33,] -0.34664775 1.08808342
[34,] -2.36756400 -0.34664775
[35,] -2.24024316 -2.36756400
[36,] -2.09737377 -2.24024316
[37,] 0.23766532 -2.09737377
[38,] 0.02021138 0.23766532
[39,] 0.32402686 0.02021138
[40,] -0.58773077 0.32402686
[41,] -0.23923214 -0.58773077
[42,] -0.68394907 -0.23923214
[43,] -0.86852523 -0.68394907
[44,] -1.05332641 -0.86852523
[45,] -1.56983698 -1.05332641
[46,] -2.02197225 -1.56983698
[47,] -2.32789171 -2.02197225
[48,] -2.56036494 -2.32789171
[49,] -2.97076728 -2.56036494
[50,] -3.20699706 -2.97076728
[51,] 2.32055295 -3.20699706
[52,] 2.99274091 2.32055295
[53,] 3.91024609 2.99274091
[54,] 4.05832184 3.91024609
[55,] 3.07375291 4.05832184
[56,] 2.15788103 3.07375291
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.62602641 0.77542748
2 0.64524269 0.62602641
3 0.24566673 0.64524269
4 1.14548149 0.24566673
5 0.61321216 1.14548149
6 0.47824484 0.61321216
7 0.43160550 0.47824484
8 -0.06530269 0.43160550
9 -0.50875669 -0.06530269
10 -0.37866407 -0.50875669
11 -0.04652853 -0.37866407
12 0.25407505 -0.04652853
13 0.99813825 0.25407505
14 0.87026645 0.99813825
15 1.41359163 0.87026645
16 -0.05210644 1.41359163
17 -0.16763258 -0.05210644
18 0.17604579 -0.16763258
19 -0.79848666 0.17604579
20 -1.97043816 -0.79848666
21 -1.29426877 -1.97043816
22 -2.71885063 -1.29426877
23 -1.77197703 -2.71885063
24 0.46238709 -1.77197703
25 0.67765345 0.46238709
26 0.97693286 0.67765345
27 0.20199815 0.97693286
28 -0.06327032 0.20199815
29 0.04685824 -0.06327032
30 2.09983719 0.04685824
31 1.65653095 2.09983719
32 1.08808342 1.65653095
33 -0.34664775 1.08808342
34 -2.36756400 -0.34664775
35 -2.24024316 -2.36756400
36 -2.09737377 -2.24024316
37 0.23766532 -2.09737377
38 0.02021138 0.23766532
39 0.32402686 0.02021138
40 -0.58773077 0.32402686
41 -0.23923214 -0.58773077
42 -0.68394907 -0.23923214
43 -0.86852523 -0.68394907
44 -1.05332641 -0.86852523
45 -1.56983698 -1.05332641
46 -2.02197225 -1.56983698
47 -2.32789171 -2.02197225
48 -2.56036494 -2.32789171
49 -2.97076728 -2.56036494
50 -3.20699706 -2.97076728
51 2.32055295 -3.20699706
52 2.99274091 2.32055295
53 3.91024609 2.99274091
54 4.05832184 3.91024609
55 3.07375291 4.05832184
56 2.15788103 3.07375291
> 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/rcomp/tmp/78z3e1292931240.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/www/rcomp/tmp/88z3e1292931240.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/www/rcomp/tmp/98z3e1292931240.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/www/rcomp/tmp/10083z1292931240.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/www/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/11491n1292931240.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/rcomp/tmp/1279ib1292931240.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/rcomp/tmp/133jx21292931240.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/rcomp/tmp/14p2wq1292931240.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/rcomp/tmp/15zbds1292931240.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/rcomp/tmp/16v3bj1292931240.tab")
+ }
> try(system("convert tmp/1cpon1292931240.ps tmp/1cpon1292931240.png",intern=TRUE))
character(0)
> try(system("convert tmp/2cpon1292931240.ps tmp/2cpon1292931240.png",intern=TRUE))
character(0)
> try(system("convert tmp/34gn81292931240.ps tmp/34gn81292931240.png",intern=TRUE))
character(0)
> try(system("convert tmp/44gn81292931240.ps tmp/44gn81292931240.png",intern=TRUE))
character(0)
> try(system("convert tmp/54gn81292931240.ps tmp/54gn81292931240.png",intern=TRUE))
character(0)
> try(system("convert tmp/6f84t1292931240.ps tmp/6f84t1292931240.png",intern=TRUE))
character(0)
> try(system("convert tmp/78z3e1292931240.ps tmp/78z3e1292931240.png",intern=TRUE))
character(0)
> try(system("convert tmp/88z3e1292931240.ps tmp/88z3e1292931240.png",intern=TRUE))
character(0)
> try(system("convert tmp/98z3e1292931240.ps tmp/98z3e1292931240.png",intern=TRUE))
character(0)
> try(system("convert tmp/10083z1292931240.ps tmp/10083z1292931240.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.090 0.720 3.807