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(106370
+ ,100.3
+ ,109375
+ ,101.9
+ ,116476
+ ,102.1
+ ,123297
+ ,103.2
+ ,114813
+ ,103.7
+ ,117925
+ ,106.2
+ ,126466
+ ,107.7
+ ,131235
+ ,109.9
+ ,120546
+ ,111.7
+ ,123791
+ ,114.9
+ ,129813
+ ,116
+ ,133463
+ ,118.3
+ ,122987
+ ,120.4
+ ,125418
+ ,126
+ ,130199
+ ,128.1
+ ,133016
+ ,130.1
+ ,121454
+ ,130.8
+ ,122044
+ ,133.6
+ ,128313
+ ,134.2
+ ,131556
+ ,135.5
+ ,120027
+ ,136.2
+ ,123001
+ ,139.1
+ ,130111
+ ,139
+ ,132524
+ ,139.6
+ ,123742
+ ,138.7
+ ,124931
+ ,140.9
+ ,133646
+ ,141.3
+ ,136557
+ ,141.8
+ ,127509
+ ,142
+ ,128945
+ ,144.5
+ ,137191
+ ,144.6
+ ,139716
+ ,145.5
+ ,129083
+ ,146.8
+ ,131604
+ ,149.5
+ ,139413
+ ,149.9
+ ,143125
+ ,150.1
+ ,133948
+ ,150.9
+ ,137116
+ ,152.8
+ ,144864
+ ,153.1
+ ,149277
+ ,154
+ ,138796
+ ,154.9
+ ,143258
+ ,156.9
+ ,150034
+ ,158.4
+ ,154708
+ ,159.7
+ ,144888
+ ,160.2
+ ,148762
+ ,163.2
+ ,156500
+ ,163.7
+ ,161088
+ ,164.4
+ ,152772
+ ,163.7
+ ,158011
+ ,165.5
+ ,163318
+ ,165.6
+ ,169969
+ ,166.8
+ ,162269
+ ,167.5
+ ,165765
+ ,170.6
+ ,170600
+ ,170.9
+ ,174681
+ ,172
+ ,166364
+ ,171.8
+ ,170240
+ ,173.9
+ ,176150
+ ,174
+ ,182056
+ ,173.8
+ ,172218
+ ,173.9
+ ,177856
+ ,176
+ ,182253
+ ,176.6
+ ,188090
+ ,178.2
+ ,176863
+ ,179.2
+ ,183273
+ ,181.3
+ ,187969
+ ,181.8
+ ,194650
+ ,182.9
+ ,183036
+ ,183.8
+ ,189516
+ ,186.3
+ ,193805
+ ,187.4
+ ,200499
+ ,189.2
+ ,188142
+ ,189.7
+ ,193732
+ ,191.9
+ ,197126
+ ,192.6
+ ,205140
+ ,193.7
+ ,191751
+ ,194.2
+ ,196700
+ ,197.6
+ ,199784
+ ,199.3
+ ,207360
+ ,201.4
+ ,196101
+ ,203
+ ,200824
+ ,206.3
+ ,205743
+ ,207.1
+ ,212489
+ ,209.8
+ ,200810
+ ,211.1
+ ,203683
+ ,215.3
+ ,207286
+ ,217.4
+ ,210910
+ ,215.5
+ ,194915
+ ,210.9
+ ,217920
+ ,212.6)
+ ,dim=c(2
+ ,90)
+ ,dimnames=list(c('HFCE'
+ ,'RPI')
+ ,1:90))
> y <- array(NA,dim=c(2,90),dimnames=list(c('HFCE','RPI'),1:90))
> 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'
> #'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
HFCE RPI
1 106370 100.3
2 109375 101.9
3 116476 102.1
4 123297 103.2
5 114813 103.7
6 117925 106.2
7 126466 107.7
8 131235 109.9
9 120546 111.7
10 123791 114.9
11 129813 116.0
12 133463 118.3
13 122987 120.4
14 125418 126.0
15 130199 128.1
16 133016 130.1
17 121454 130.8
18 122044 133.6
19 128313 134.2
20 131556 135.5
21 120027 136.2
22 123001 139.1
23 130111 139.0
24 132524 139.6
25 123742 138.7
26 124931 140.9
27 133646 141.3
28 136557 141.8
29 127509 142.0
30 128945 144.5
31 137191 144.6
32 139716 145.5
33 129083 146.8
34 131604 149.5
35 139413 149.9
36 143125 150.1
37 133948 150.9
38 137116 152.8
39 144864 153.1
40 149277 154.0
41 138796 154.9
42 143258 156.9
43 150034 158.4
44 154708 159.7
45 144888 160.2
46 148762 163.2
47 156500 163.7
48 161088 164.4
49 152772 163.7
50 158011 165.5
51 163318 165.6
52 169969 166.8
53 162269 167.5
54 165765 170.6
55 170600 170.9
56 174681 172.0
57 166364 171.8
58 170240 173.9
59 176150 174.0
60 182056 173.8
61 172218 173.9
62 177856 176.0
63 182253 176.6
64 188090 178.2
65 176863 179.2
66 183273 181.3
67 187969 181.8
68 194650 182.9
69 183036 183.8
70 189516 186.3
71 193805 187.4
72 200499 189.2
73 188142 189.7
74 193732 191.9
75 197126 192.6
76 205140 193.7
77 191751 194.2
78 196700 197.6
79 199784 199.3
80 207360 201.4
81 196101 203.0
82 200824 206.3
83 205743 207.1
84 212489 209.8
85 200810 211.1
86 203683 215.3
87 207286 217.4
88 210910 215.5
89 194915 210.9
90 217920 212.6
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) RPI
7117.6 942.1
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-16360.1 -6304.3 -110.2 6725.4 20578.8
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7117.58 5207.51 1.367 0.175
RPI 942.12 31.94 29.494 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9504 on 88 degrees of freedom
Multiple R-squared: 0.9081, Adjusted R-squared: 0.9071
F-statistic: 869.9 on 1 and 88 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.16856027 0.33712055 0.831439726
[2,] 0.14133657 0.28267315 0.858663426
[3,] 0.08060798 0.16121597 0.919392017
[4,] 0.04885611 0.09771221 0.951143895
[5,] 0.14385497 0.28770993 0.856145035
[6,] 0.15495594 0.30991188 0.845044062
[7,] 0.12620784 0.25241568 0.873792160
[8,] 0.11938495 0.23876989 0.880615053
[9,] 0.22179237 0.44358474 0.778207630
[10,] 0.27616878 0.55233757 0.723831216
[11,] 0.24291860 0.48583720 0.757081399
[12,] 0.21135224 0.42270449 0.788647756
[13,] 0.31676051 0.63352101 0.683239495
[14,] 0.36235802 0.72471604 0.637641981
[15,] 0.30048946 0.60097892 0.699510538
[16,] 0.24734227 0.49468453 0.752657733
[17,] 0.31098532 0.62197065 0.689014677
[18,] 0.30631951 0.61263901 0.693680494
[19,] 0.24573406 0.49146813 0.754265936
[20,] 0.20145847 0.40291694 0.798541528
[21,] 0.18379016 0.36758032 0.816209842
[22,] 0.16293931 0.32587862 0.837060691
[23,] 0.13592733 0.27185465 0.864072673
[24,] 0.12947444 0.25894888 0.870525558
[25,] 0.10529006 0.21058011 0.894709944
[26,] 0.08609179 0.17218358 0.913908209
[27,] 0.07779656 0.15559313 0.922203436
[28,] 0.07968750 0.15937500 0.920312502
[29,] 0.07448375 0.14896750 0.925516249
[30,] 0.07019639 0.14039278 0.929803609
[31,] 0.06431607 0.12863213 0.935683933
[32,] 0.07142684 0.14285369 0.928573157
[33,] 0.06934675 0.13869350 0.930653249
[34,] 0.06801524 0.13603048 0.931984761
[35,] 0.07693074 0.15386149 0.923069257
[36,] 0.10906658 0.21813317 0.890933415
[37,] 0.11948120 0.23896240 0.880518798
[38,] 0.13192817 0.26385634 0.868071831
[39,] 0.16055481 0.32110963 0.839445186
[40,] 0.21880171 0.43760341 0.781198295
[41,] 0.27839833 0.55679666 0.721601671
[42,] 0.37235554 0.74471108 0.627644461
[43,] 0.45844882 0.91689764 0.541551179
[44,] 0.56174337 0.87651326 0.438256628
[45,] 0.65867050 0.68265899 0.341329496
[46,] 0.73956247 0.52087506 0.260437530
[47,] 0.80363914 0.39272173 0.196360864
[48,] 0.87563678 0.24872643 0.124363215
[49,] 0.90532870 0.18934260 0.094671298
[50,] 0.93158552 0.13682896 0.068414482
[51,] 0.94651240 0.10697521 0.053487604
[52,] 0.95802737 0.08394527 0.041972633
[53,] 0.97342002 0.05315997 0.026579983
[54,] 0.98250539 0.03498922 0.017494611
[55,] 0.98500898 0.02998204 0.014991022
[56,] 0.98871260 0.02257481 0.011287403
[57,] 0.99170874 0.01658252 0.008291262
[58,] 0.99210723 0.01578553 0.007892767
[59,] 0.99160837 0.01678326 0.008391630
[60,] 0.99244770 0.01510460 0.007552301
[61,] 0.99455791 0.01088417 0.005442087
[62,] 0.99394221 0.01211557 0.006057786
[63,] 0.99217084 0.01565832 0.007829158
[64,] 0.99281546 0.01436909 0.007184543
[65,] 0.99286286 0.01427428 0.007137141
[66,] 0.98986340 0.02027321 0.010136603
[67,] 0.98497105 0.03005790 0.015028950
[68,] 0.98580369 0.02839262 0.014196309
[69,] 0.98248795 0.03502410 0.017512051
[70,] 0.97256229 0.05487543 0.027437714
[71,] 0.95640962 0.08718077 0.043590383
[72,] 0.96556650 0.06886701 0.034433503
[73,] 0.94884531 0.10230937 0.051154686
[74,] 0.91817716 0.16364569 0.081822845
[75,] 0.86960322 0.26079355 0.130396777
[76,] 0.85787592 0.28424816 0.142124082
[77,] 0.79911529 0.40176941 0.200884707
[78,] 0.71073339 0.57853322 0.289266612
[79,] 0.58474844 0.83050312 0.415251561
[80,] 0.56824190 0.86351621 0.431758104
[81,] 0.41062877 0.82125754 0.589371231
> postscript(file="/var/www/html/rcomp/tmp/1ed1u1258722754.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/2bowv1258722754.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/3u0bg1258722754.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/4onp71258722754.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/5s8yo1258722754.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 = 90
Frequency = 1
1 2 3 4 5 6
4758.0824 6255.6951 13168.2717 18952.9430 9997.8845 10754.5919
7 8 9 10 11 12
17882.4163 20578.7589 8193.9482 8424.1737 13409.8450 14892.9758
13 14 15 16 17 18
2438.5301 -406.3253 2396.2289 3328.9949 -8892.4870 -10940.4147
19 20 21 22 23 24
-5236.6849 -3218.4371 -15406.9190 -15165.0584 -7960.8467 -6113.1169
25 26 27 28 29 30
-14047.2116 -14930.8690 -6592.7158 -4152.7744 -13389.1978 -14308.4903
31 32 33 34 35 36
-6156.7020 -4479.6074 -16337.3595 -16360.0755 -8927.9223 -5404.3457
37 38 39 40 41 42
-15335.0393 -13957.0617 -6491.6968 -2926.6021 -14255.5075 -11677.7415
43 44 45 46 47 48
-6314.9171 -2865.6692 -13156.7277 -12109.0788 -4842.1373 -913.6192
49 50 51 52 53 54
-8570.1373 -5026.9480 185.8403 5706.2999 -2653.1820 -2077.7448
55 56 57 58 59 60
2474.6201 5519.2913 -2609.2853 -711.7310 5104.0573 11198.4807
61 62 63 64 65 66
1266.2690 4925.8232 8757.5530 13087.1658 918.0487 5349.6030
67 68 69 70 71 72
9574.5445 15219.2157 2757.3104 6882.0178 10134.6891 15132.8784
73 74 75 76 77 78
2304.8199 5822.1625 8556.6805 15534.3518 1674.2933 3420.0954
79 80 81 82 83 84
4902.4964 10500.0507 -2266.3366 -652.3228 3512.9836 7715.2676
85 86 87 88 89 90
-5188.4845 -6272.3760 -4647.8218 766.2006 -10895.0611 10508.3399
> postscript(file="/var/www/html/rcomp/tmp/6wy071258722754.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 = 90
Frequency = 1
lag(myerror, k = 1) myerror
0 4758.0824 NA
1 6255.6951 4758.0824
2 13168.2717 6255.6951
3 18952.9430 13168.2717
4 9997.8845 18952.9430
5 10754.5919 9997.8845
6 17882.4163 10754.5919
7 20578.7589 17882.4163
8 8193.9482 20578.7589
9 8424.1737 8193.9482
10 13409.8450 8424.1737
11 14892.9758 13409.8450
12 2438.5301 14892.9758
13 -406.3253 2438.5301
14 2396.2289 -406.3253
15 3328.9949 2396.2289
16 -8892.4870 3328.9949
17 -10940.4147 -8892.4870
18 -5236.6849 -10940.4147
19 -3218.4371 -5236.6849
20 -15406.9190 -3218.4371
21 -15165.0584 -15406.9190
22 -7960.8467 -15165.0584
23 -6113.1169 -7960.8467
24 -14047.2116 -6113.1169
25 -14930.8690 -14047.2116
26 -6592.7158 -14930.8690
27 -4152.7744 -6592.7158
28 -13389.1978 -4152.7744
29 -14308.4903 -13389.1978
30 -6156.7020 -14308.4903
31 -4479.6074 -6156.7020
32 -16337.3595 -4479.6074
33 -16360.0755 -16337.3595
34 -8927.9223 -16360.0755
35 -5404.3457 -8927.9223
36 -15335.0393 -5404.3457
37 -13957.0617 -15335.0393
38 -6491.6968 -13957.0617
39 -2926.6021 -6491.6968
40 -14255.5075 -2926.6021
41 -11677.7415 -14255.5075
42 -6314.9171 -11677.7415
43 -2865.6692 -6314.9171
44 -13156.7277 -2865.6692
45 -12109.0788 -13156.7277
46 -4842.1373 -12109.0788
47 -913.6192 -4842.1373
48 -8570.1373 -913.6192
49 -5026.9480 -8570.1373
50 185.8403 -5026.9480
51 5706.2999 185.8403
52 -2653.1820 5706.2999
53 -2077.7448 -2653.1820
54 2474.6201 -2077.7448
55 5519.2913 2474.6201
56 -2609.2853 5519.2913
57 -711.7310 -2609.2853
58 5104.0573 -711.7310
59 11198.4807 5104.0573
60 1266.2690 11198.4807
61 4925.8232 1266.2690
62 8757.5530 4925.8232
63 13087.1658 8757.5530
64 918.0487 13087.1658
65 5349.6030 918.0487
66 9574.5445 5349.6030
67 15219.2157 9574.5445
68 2757.3104 15219.2157
69 6882.0178 2757.3104
70 10134.6891 6882.0178
71 15132.8784 10134.6891
72 2304.8199 15132.8784
73 5822.1625 2304.8199
74 8556.6805 5822.1625
75 15534.3518 8556.6805
76 1674.2933 15534.3518
77 3420.0954 1674.2933
78 4902.4964 3420.0954
79 10500.0507 4902.4964
80 -2266.3366 10500.0507
81 -652.3228 -2266.3366
82 3512.9836 -652.3228
83 7715.2676 3512.9836
84 -5188.4845 7715.2676
85 -6272.3760 -5188.4845
86 -4647.8218 -6272.3760
87 766.2006 -4647.8218
88 -10895.0611 766.2006
89 10508.3399 -10895.0611
90 NA 10508.3399
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 6255.6951 4758.0824
[2,] 13168.2717 6255.6951
[3,] 18952.9430 13168.2717
[4,] 9997.8845 18952.9430
[5,] 10754.5919 9997.8845
[6,] 17882.4163 10754.5919
[7,] 20578.7589 17882.4163
[8,] 8193.9482 20578.7589
[9,] 8424.1737 8193.9482
[10,] 13409.8450 8424.1737
[11,] 14892.9758 13409.8450
[12,] 2438.5301 14892.9758
[13,] -406.3253 2438.5301
[14,] 2396.2289 -406.3253
[15,] 3328.9949 2396.2289
[16,] -8892.4870 3328.9949
[17,] -10940.4147 -8892.4870
[18,] -5236.6849 -10940.4147
[19,] -3218.4371 -5236.6849
[20,] -15406.9190 -3218.4371
[21,] -15165.0584 -15406.9190
[22,] -7960.8467 -15165.0584
[23,] -6113.1169 -7960.8467
[24,] -14047.2116 -6113.1169
[25,] -14930.8690 -14047.2116
[26,] -6592.7158 -14930.8690
[27,] -4152.7744 -6592.7158
[28,] -13389.1978 -4152.7744
[29,] -14308.4903 -13389.1978
[30,] -6156.7020 -14308.4903
[31,] -4479.6074 -6156.7020
[32,] -16337.3595 -4479.6074
[33,] -16360.0755 -16337.3595
[34,] -8927.9223 -16360.0755
[35,] -5404.3457 -8927.9223
[36,] -15335.0393 -5404.3457
[37,] -13957.0617 -15335.0393
[38,] -6491.6968 -13957.0617
[39,] -2926.6021 -6491.6968
[40,] -14255.5075 -2926.6021
[41,] -11677.7415 -14255.5075
[42,] -6314.9171 -11677.7415
[43,] -2865.6692 -6314.9171
[44,] -13156.7277 -2865.6692
[45,] -12109.0788 -13156.7277
[46,] -4842.1373 -12109.0788
[47,] -913.6192 -4842.1373
[48,] -8570.1373 -913.6192
[49,] -5026.9480 -8570.1373
[50,] 185.8403 -5026.9480
[51,] 5706.2999 185.8403
[52,] -2653.1820 5706.2999
[53,] -2077.7448 -2653.1820
[54,] 2474.6201 -2077.7448
[55,] 5519.2913 2474.6201
[56,] -2609.2853 5519.2913
[57,] -711.7310 -2609.2853
[58,] 5104.0573 -711.7310
[59,] 11198.4807 5104.0573
[60,] 1266.2690 11198.4807
[61,] 4925.8232 1266.2690
[62,] 8757.5530 4925.8232
[63,] 13087.1658 8757.5530
[64,] 918.0487 13087.1658
[65,] 5349.6030 918.0487
[66,] 9574.5445 5349.6030
[67,] 15219.2157 9574.5445
[68,] 2757.3104 15219.2157
[69,] 6882.0178 2757.3104
[70,] 10134.6891 6882.0178
[71,] 15132.8784 10134.6891
[72,] 2304.8199 15132.8784
[73,] 5822.1625 2304.8199
[74,] 8556.6805 5822.1625
[75,] 15534.3518 8556.6805
[76,] 1674.2933 15534.3518
[77,] 3420.0954 1674.2933
[78,] 4902.4964 3420.0954
[79,] 10500.0507 4902.4964
[80,] -2266.3366 10500.0507
[81,] -652.3228 -2266.3366
[82,] 3512.9836 -652.3228
[83,] 7715.2676 3512.9836
[84,] -5188.4845 7715.2676
[85,] -6272.3760 -5188.4845
[86,] -4647.8218 -6272.3760
[87,] 766.2006 -4647.8218
[88,] -10895.0611 766.2006
[89,] 10508.3399 -10895.0611
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 6255.6951 4758.0824
2 13168.2717 6255.6951
3 18952.9430 13168.2717
4 9997.8845 18952.9430
5 10754.5919 9997.8845
6 17882.4163 10754.5919
7 20578.7589 17882.4163
8 8193.9482 20578.7589
9 8424.1737 8193.9482
10 13409.8450 8424.1737
11 14892.9758 13409.8450
12 2438.5301 14892.9758
13 -406.3253 2438.5301
14 2396.2289 -406.3253
15 3328.9949 2396.2289
16 -8892.4870 3328.9949
17 -10940.4147 -8892.4870
18 -5236.6849 -10940.4147
19 -3218.4371 -5236.6849
20 -15406.9190 -3218.4371
21 -15165.0584 -15406.9190
22 -7960.8467 -15165.0584
23 -6113.1169 -7960.8467
24 -14047.2116 -6113.1169
25 -14930.8690 -14047.2116
26 -6592.7158 -14930.8690
27 -4152.7744 -6592.7158
28 -13389.1978 -4152.7744
29 -14308.4903 -13389.1978
30 -6156.7020 -14308.4903
31 -4479.6074 -6156.7020
32 -16337.3595 -4479.6074
33 -16360.0755 -16337.3595
34 -8927.9223 -16360.0755
35 -5404.3457 -8927.9223
36 -15335.0393 -5404.3457
37 -13957.0617 -15335.0393
38 -6491.6968 -13957.0617
39 -2926.6021 -6491.6968
40 -14255.5075 -2926.6021
41 -11677.7415 -14255.5075
42 -6314.9171 -11677.7415
43 -2865.6692 -6314.9171
44 -13156.7277 -2865.6692
45 -12109.0788 -13156.7277
46 -4842.1373 -12109.0788
47 -913.6192 -4842.1373
48 -8570.1373 -913.6192
49 -5026.9480 -8570.1373
50 185.8403 -5026.9480
51 5706.2999 185.8403
52 -2653.1820 5706.2999
53 -2077.7448 -2653.1820
54 2474.6201 -2077.7448
55 5519.2913 2474.6201
56 -2609.2853 5519.2913
57 -711.7310 -2609.2853
58 5104.0573 -711.7310
59 11198.4807 5104.0573
60 1266.2690 11198.4807
61 4925.8232 1266.2690
62 8757.5530 4925.8232
63 13087.1658 8757.5530
64 918.0487 13087.1658
65 5349.6030 918.0487
66 9574.5445 5349.6030
67 15219.2157 9574.5445
68 2757.3104 15219.2157
69 6882.0178 2757.3104
70 10134.6891 6882.0178
71 15132.8784 10134.6891
72 2304.8199 15132.8784
73 5822.1625 2304.8199
74 8556.6805 5822.1625
75 15534.3518 8556.6805
76 1674.2933 15534.3518
77 3420.0954 1674.2933
78 4902.4964 3420.0954
79 10500.0507 4902.4964
80 -2266.3366 10500.0507
81 -652.3228 -2266.3366
82 3512.9836 -652.3228
83 7715.2676 3512.9836
84 -5188.4845 7715.2676
85 -6272.3760 -5188.4845
86 -4647.8218 -6272.3760
87 766.2006 -4647.8218
88 -10895.0611 766.2006
89 10508.3399 -10895.0611
> 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/7poba1258722754.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/8ca0n1258722754.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/9yexz1258722754.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/10pt7z1258722754.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/11m9mc1258722754.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/12lalj1258722754.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/13tezc1258722754.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/14sar61258722754.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/15wnwv1258722754.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/160yhd1258722754.tab")
+ }
>
> system("convert tmp/1ed1u1258722754.ps tmp/1ed1u1258722754.png")
> system("convert tmp/2bowv1258722754.ps tmp/2bowv1258722754.png")
> system("convert tmp/3u0bg1258722754.ps tmp/3u0bg1258722754.png")
> system("convert tmp/4onp71258722754.ps tmp/4onp71258722754.png")
> system("convert tmp/5s8yo1258722754.ps tmp/5s8yo1258722754.png")
> system("convert tmp/6wy071258722754.ps tmp/6wy071258722754.png")
> system("convert tmp/7poba1258722754.ps tmp/7poba1258722754.png")
> system("convert tmp/8ca0n1258722754.ps tmp/8ca0n1258722754.png")
> system("convert tmp/9yexz1258722754.ps tmp/9yexz1258722754.png")
> system("convert tmp/10pt7z1258722754.ps tmp/10pt7z1258722754.png")
>
>
> proc.time()
user system elapsed
2.703 1.555 3.119