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(519164,0.9,517009,1.3,509933,1.4,509127,1.5,500857,1.1,506971,1.6,569323,1.5,579714,1.6,577992,1.7,565464,1.6,547344,1.7,554788,1.6,562325,1.6,560854,1.3,555332,1.1,543599,1.6,536662,1.9,542722,1.6,593530,1.7,610763,1.6,612613,1.4,611324,2.1,594167,1.9,595454,1.7,590865,1.8,589379,2,584428,2.5,573100,2.1,567456,2.1,569028,2.3,620735,2.4,628884,2.4,628232,2.3,612117,1.7,595404,2,597141,2.3,593408,2,590072,2,579799,1.3,574205,1.7,572775,1.9,572942,1.7,619567,1.6,625809,1.7,619916,1.8,587625,1.9,565742,1.9,557274,1.9,560576,2,548854,2.1,531673,1.9,525919,1.9,511038,1.3,498662,1.3,555362,1.4,564591,1.2,541657,1.3,527070,1.8,509846,2.2,514258,2.6,516922,2.8,507561,3.1,492622,3.9,490243,3.7,469357,4.6,477580,5.1,528379,5.2,533590,4.9,517945,5.1,506174,4.8,501866,3.9,516141,3.5),dim=c(2,72),dimnames=list(c('TWIB','GI'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('TWIB','GI'),1:72))
> 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
TWIB GI M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 519164 0.9 1 0 0 0 0 0 0 0 0 0 0 1
2 517009 1.3 0 1 0 0 0 0 0 0 0 0 0 2
3 509933 1.4 0 0 1 0 0 0 0 0 0 0 0 3
4 509127 1.5 0 0 0 1 0 0 0 0 0 0 0 4
5 500857 1.1 0 0 0 0 1 0 0 0 0 0 0 5
6 506971 1.6 0 0 0 0 0 1 0 0 0 0 0 6
7 569323 1.5 0 0 0 0 0 0 1 0 0 0 0 7
8 579714 1.6 0 0 0 0 0 0 0 1 0 0 0 8
9 577992 1.7 0 0 0 0 0 0 0 0 1 0 0 9
10 565464 1.6 0 0 0 0 0 0 0 0 0 1 0 10
11 547344 1.7 0 0 0 0 0 0 0 0 0 0 1 11
12 554788 1.6 0 0 0 0 0 0 0 0 0 0 0 12
13 562325 1.6 1 0 0 0 0 0 0 0 0 0 0 13
14 560854 1.3 0 1 0 0 0 0 0 0 0 0 0 14
15 555332 1.1 0 0 1 0 0 0 0 0 0 0 0 15
16 543599 1.6 0 0 0 1 0 0 0 0 0 0 0 16
17 536662 1.9 0 0 0 0 1 0 0 0 0 0 0 17
18 542722 1.6 0 0 0 0 0 1 0 0 0 0 0 18
19 593530 1.7 0 0 0 0 0 0 1 0 0 0 0 19
20 610763 1.6 0 0 0 0 0 0 0 1 0 0 0 20
21 612613 1.4 0 0 0 0 0 0 0 0 1 0 0 21
22 611324 2.1 0 0 0 0 0 0 0 0 0 1 0 22
23 594167 1.9 0 0 0 0 0 0 0 0 0 0 1 23
24 595454 1.7 0 0 0 0 0 0 0 0 0 0 0 24
25 590865 1.8 1 0 0 0 0 0 0 0 0 0 0 25
26 589379 2.0 0 1 0 0 0 0 0 0 0 0 0 26
27 584428 2.5 0 0 1 0 0 0 0 0 0 0 0 27
28 573100 2.1 0 0 0 1 0 0 0 0 0 0 0 28
29 567456 2.1 0 0 0 0 1 0 0 0 0 0 0 29
30 569028 2.3 0 0 0 0 0 1 0 0 0 0 0 30
31 620735 2.4 0 0 0 0 0 0 1 0 0 0 0 31
32 628884 2.4 0 0 0 0 0 0 0 1 0 0 0 32
33 628232 2.3 0 0 0 0 0 0 0 0 1 0 0 33
34 612117 1.7 0 0 0 0 0 0 0 0 0 1 0 34
35 595404 2.0 0 0 0 0 0 0 0 0 0 0 1 35
36 597141 2.3 0 0 0 0 0 0 0 0 0 0 0 36
37 593408 2.0 1 0 0 0 0 0 0 0 0 0 0 37
38 590072 2.0 0 1 0 0 0 0 0 0 0 0 0 38
39 579799 1.3 0 0 1 0 0 0 0 0 0 0 0 39
40 574205 1.7 0 0 0 1 0 0 0 0 0 0 0 40
41 572775 1.9 0 0 0 0 1 0 0 0 0 0 0 41
42 572942 1.7 0 0 0 0 0 1 0 0 0 0 0 42
43 619567 1.6 0 0 0 0 0 0 1 0 0 0 0 43
44 625809 1.7 0 0 0 0 0 0 0 1 0 0 0 44
45 619916 1.8 0 0 0 0 0 0 0 0 1 0 0 45
46 587625 1.9 0 0 0 0 0 0 0 0 0 1 0 46
47 565742 1.9 0 0 0 0 0 0 0 0 0 0 1 47
48 557274 1.9 0 0 0 0 0 0 0 0 0 0 0 48
49 560576 2.0 1 0 0 0 0 0 0 0 0 0 0 49
50 548854 2.1 0 1 0 0 0 0 0 0 0 0 0 50
51 531673 1.9 0 0 1 0 0 0 0 0 0 0 0 51
52 525919 1.9 0 0 0 1 0 0 0 0 0 0 0 52
53 511038 1.3 0 0 0 0 1 0 0 0 0 0 0 53
54 498662 1.3 0 0 0 0 0 1 0 0 0 0 0 54
55 555362 1.4 0 0 0 0 0 0 1 0 0 0 0 55
56 564591 1.2 0 0 0 0 0 0 0 1 0 0 0 56
57 541657 1.3 0 0 0 0 0 0 0 0 1 0 0 57
58 527070 1.8 0 0 0 0 0 0 0 0 0 1 0 58
59 509846 2.2 0 0 0 0 0 0 0 0 0 0 1 59
60 514258 2.6 0 0 0 0 0 0 0 0 0 0 0 60
61 516922 2.8 1 0 0 0 0 0 0 0 0 0 0 61
62 507561 3.1 0 1 0 0 0 0 0 0 0 0 0 62
63 492622 3.9 0 0 1 0 0 0 0 0 0 0 0 63
64 490243 3.7 0 0 0 1 0 0 0 0 0 0 0 64
65 469357 4.6 0 0 0 0 1 0 0 0 0 0 0 65
66 477580 5.1 0 0 0 0 0 1 0 0 0 0 0 66
67 528379 5.2 0 0 0 0 0 0 1 0 0 0 0 67
68 533590 4.9 0 0 0 0 0 0 0 1 0 0 0 68
69 517945 5.1 0 0 0 0 0 0 0 0 1 0 0 69
70 506174 4.8 0 0 0 0 0 0 0 0 0 1 0 70
71 501866 3.9 0 0 0 0 0 0 0 0 0 0 1 71
72 516141 3.5 0 0 0 0 0 0 0 0 0 0 0 72
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) GI M1 M2 M3 M4
598405.6 -13369.7 -7413.9 -10484.1 -19514.0 -24596.5
M5 M6 M7 M8 M9 M10
-33088.0 -29609.7 24293.0 33102.7 26340.9 12537.8
M11 t
-3739.7 -291.9
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-59503 -21223 -3958 31804 46841
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 598405.6 16560.0 36.136 <2e-16 ***
GI -13369.7 5163.8 -2.589 0.0121 *
M1 -7413.9 19487.0 -0.380 0.7050
M2 -10484.1 19465.4 -0.539 0.5922
M3 -19514.0 19448.1 -1.003 0.3198
M4 -24596.5 19434.7 -1.266 0.2107
M5 -33088.0 19424.8 -1.703 0.0938 .
M6 -29609.7 19430.2 -1.524 0.1330
M7 24293.0 19419.8 1.251 0.2160
M8 33102.7 19391.0 1.707 0.0931 .
M9 26340.9 19384.3 1.359 0.1794
M10 12537.8 19382.0 0.647 0.5203
M11 -3739.7 19370.7 -0.193 0.8476
t -291.9 256.8 -1.136 0.2605
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 33550 on 58 degrees of freedom
Multiple R-squared: 0.451, Adjusted R-squared: 0.328
F-statistic: 3.665 on 13 and 58 DF, p-value: 0.0003044
> 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,] 8.087245e-03 0.0161744891 9.919128e-01
[2,] 2.656475e-03 0.0053129498 9.973435e-01
[3,] 8.036806e-03 0.0160736122 9.919632e-01
[4,] 4.656379e-03 0.0093127579 9.953436e-01
[5,] 2.011156e-03 0.0040223111 9.979888e-01
[6,] 1.253954e-03 0.0025079085 9.987460e-01
[7,] 8.896572e-04 0.0017793144 9.991103e-01
[8,] 4.172829e-04 0.0008345658 9.995827e-01
[9,] 3.514721e-04 0.0007029443 9.996485e-01
[10,] 2.150572e-04 0.0004301144 9.997849e-01
[11,] 8.631594e-05 0.0001726319 9.999137e-01
[12,] 8.793618e-05 0.0001758724 9.999121e-01
[13,] 5.571143e-05 0.0001114229 9.999443e-01
[14,] 5.957572e-05 0.0001191514 9.999404e-01
[15,] 1.368180e-04 0.0002736360 9.998632e-01
[16,] 7.413554e-04 0.0014827109 9.992586e-01
[17,] 1.594493e-03 0.0031889861 9.984055e-01
[18,] 1.399287e-02 0.0279857361 9.860071e-01
[19,] 3.115379e-02 0.0623075713 9.688462e-01
[20,] 8.482063e-02 0.1696412689 9.151794e-01
[21,] 1.353402e-01 0.2706803496 8.646598e-01
[22,] 1.617735e-01 0.3235470983 8.382265e-01
[23,] 1.479022e-01 0.2958043010 8.520978e-01
[24,] 1.220074e-01 0.2440148086 8.779926e-01
[25,] 1.048196e-01 0.2096391532 8.951804e-01
[26,] 1.120285e-01 0.2240570655 8.879715e-01
[27,] 1.174516e-01 0.2349032112 8.825484e-01
[28,] 1.389523e-01 0.2779045698 8.610477e-01
[29,] 4.448871e-01 0.8897742782 5.551129e-01
[30,] 7.815475e-01 0.4369049976 2.184525e-01
[31,] 8.922708e-01 0.2154584178 1.077292e-01
[32,] 9.249912e-01 0.1500176262 7.500881e-02
[33,] 9.557665e-01 0.0884670148 4.423351e-02
[34,] 9.815081e-01 0.0369838054 1.849190e-02
[35,] 9.914348e-01 0.0171303272 8.565164e-03
[36,] 9.981174e-01 0.0037651687 1.882584e-03
[37,] 9.999059e-01 0.0001882995 9.414976e-05
[38,] 9.998120e-01 0.0003759752 1.879876e-04
[39,] 9.981682e-01 0.0036636931 1.831847e-03
> postscript(file="/var/www/html/rcomp/tmp/1177e1258743587.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/2p8tx1258743587.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/3fzi11258743587.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/41ja81258743587.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/5n7fg1258743587.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 = 72
Frequency = 1
1 2 3 4 5 6
-59503.1489 -52948.2459 -49365.4297 -43460.1078 -48294.6147 -38682.2459
7 8 9 10 11 12
-31278.0336 -28067.9240 -21399.2801 -21169.2288 -21382.9468 -18723.7459
13 14 15 16 17 18
-3480.9989 -5600.8563 -4474.9374 -4148.7524 1708.5011 571.1437
19 20 21 22 23 24
-894.7124 6483.4657 12713.2122 34877.9897 31616.3744 26781.6095
25 26 27 28 29 30
31235.3223 35785.2938 46840.9731 35539.4661 38678.8223 39738.2938
31 32 33 34 35 36
39171.4376 38802.5815 43867.2938 33825.5163 37692.7298 39992.7938
37 38 39 40 41 42
39954.6435 39980.6834 29670.7735 34798.9927 44826.2804 39132.8888
43 44 45 46 47 48
30810.1011 29871.2107 32368.8546 15509.8375 10196.1537 -1719.6797
49 50 51 52 53 54
10625.0331 3602.0388 -6931.0423 -7310.6861 -21430.1246 -36992.5846
55 56 57 58 59 60
-32566.4408 -34529.2285 -49072.5846 -42879.7387 -38186.5594 -31874.5296
61 62 63 64 65 66
-18830.8511 -20818.9138 -15740.3372 -15418.9126 -15488.8645 -3767.4958
67 68 69 70 71 72
-5242.3519 -12560.1054 -18477.4958 -20164.3760 -19935.7517 -14456.4481
> postscript(file="/var/www/html/rcomp/tmp/6h3o81258743587.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 -59503.1489 NA
1 -52948.2459 -59503.1489
2 -49365.4297 -52948.2459
3 -43460.1078 -49365.4297
4 -48294.6147 -43460.1078
5 -38682.2459 -48294.6147
6 -31278.0336 -38682.2459
7 -28067.9240 -31278.0336
8 -21399.2801 -28067.9240
9 -21169.2288 -21399.2801
10 -21382.9468 -21169.2288
11 -18723.7459 -21382.9468
12 -3480.9989 -18723.7459
13 -5600.8563 -3480.9989
14 -4474.9374 -5600.8563
15 -4148.7524 -4474.9374
16 1708.5011 -4148.7524
17 571.1437 1708.5011
18 -894.7124 571.1437
19 6483.4657 -894.7124
20 12713.2122 6483.4657
21 34877.9897 12713.2122
22 31616.3744 34877.9897
23 26781.6095 31616.3744
24 31235.3223 26781.6095
25 35785.2938 31235.3223
26 46840.9731 35785.2938
27 35539.4661 46840.9731
28 38678.8223 35539.4661
29 39738.2938 38678.8223
30 39171.4376 39738.2938
31 38802.5815 39171.4376
32 43867.2938 38802.5815
33 33825.5163 43867.2938
34 37692.7298 33825.5163
35 39992.7938 37692.7298
36 39954.6435 39992.7938
37 39980.6834 39954.6435
38 29670.7735 39980.6834
39 34798.9927 29670.7735
40 44826.2804 34798.9927
41 39132.8888 44826.2804
42 30810.1011 39132.8888
43 29871.2107 30810.1011
44 32368.8546 29871.2107
45 15509.8375 32368.8546
46 10196.1537 15509.8375
47 -1719.6797 10196.1537
48 10625.0331 -1719.6797
49 3602.0388 10625.0331
50 -6931.0423 3602.0388
51 -7310.6861 -6931.0423
52 -21430.1246 -7310.6861
53 -36992.5846 -21430.1246
54 -32566.4408 -36992.5846
55 -34529.2285 -32566.4408
56 -49072.5846 -34529.2285
57 -42879.7387 -49072.5846
58 -38186.5594 -42879.7387
59 -31874.5296 -38186.5594
60 -18830.8511 -31874.5296
61 -20818.9138 -18830.8511
62 -15740.3372 -20818.9138
63 -15418.9126 -15740.3372
64 -15488.8645 -15418.9126
65 -3767.4958 -15488.8645
66 -5242.3519 -3767.4958
67 -12560.1054 -5242.3519
68 -18477.4958 -12560.1054
69 -20164.3760 -18477.4958
70 -19935.7517 -20164.3760
71 -14456.4481 -19935.7517
72 NA -14456.4481
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -52948.2459 -59503.1489
[2,] -49365.4297 -52948.2459
[3,] -43460.1078 -49365.4297
[4,] -48294.6147 -43460.1078
[5,] -38682.2459 -48294.6147
[6,] -31278.0336 -38682.2459
[7,] -28067.9240 -31278.0336
[8,] -21399.2801 -28067.9240
[9,] -21169.2288 -21399.2801
[10,] -21382.9468 -21169.2288
[11,] -18723.7459 -21382.9468
[12,] -3480.9989 -18723.7459
[13,] -5600.8563 -3480.9989
[14,] -4474.9374 -5600.8563
[15,] -4148.7524 -4474.9374
[16,] 1708.5011 -4148.7524
[17,] 571.1437 1708.5011
[18,] -894.7124 571.1437
[19,] 6483.4657 -894.7124
[20,] 12713.2122 6483.4657
[21,] 34877.9897 12713.2122
[22,] 31616.3744 34877.9897
[23,] 26781.6095 31616.3744
[24,] 31235.3223 26781.6095
[25,] 35785.2938 31235.3223
[26,] 46840.9731 35785.2938
[27,] 35539.4661 46840.9731
[28,] 38678.8223 35539.4661
[29,] 39738.2938 38678.8223
[30,] 39171.4376 39738.2938
[31,] 38802.5815 39171.4376
[32,] 43867.2938 38802.5815
[33,] 33825.5163 43867.2938
[34,] 37692.7298 33825.5163
[35,] 39992.7938 37692.7298
[36,] 39954.6435 39992.7938
[37,] 39980.6834 39954.6435
[38,] 29670.7735 39980.6834
[39,] 34798.9927 29670.7735
[40,] 44826.2804 34798.9927
[41,] 39132.8888 44826.2804
[42,] 30810.1011 39132.8888
[43,] 29871.2107 30810.1011
[44,] 32368.8546 29871.2107
[45,] 15509.8375 32368.8546
[46,] 10196.1537 15509.8375
[47,] -1719.6797 10196.1537
[48,] 10625.0331 -1719.6797
[49,] 3602.0388 10625.0331
[50,] -6931.0423 3602.0388
[51,] -7310.6861 -6931.0423
[52,] -21430.1246 -7310.6861
[53,] -36992.5846 -21430.1246
[54,] -32566.4408 -36992.5846
[55,] -34529.2285 -32566.4408
[56,] -49072.5846 -34529.2285
[57,] -42879.7387 -49072.5846
[58,] -38186.5594 -42879.7387
[59,] -31874.5296 -38186.5594
[60,] -18830.8511 -31874.5296
[61,] -20818.9138 -18830.8511
[62,] -15740.3372 -20818.9138
[63,] -15418.9126 -15740.3372
[64,] -15488.8645 -15418.9126
[65,] -3767.4958 -15488.8645
[66,] -5242.3519 -3767.4958
[67,] -12560.1054 -5242.3519
[68,] -18477.4958 -12560.1054
[69,] -20164.3760 -18477.4958
[70,] -19935.7517 -20164.3760
[71,] -14456.4481 -19935.7517
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -52948.2459 -59503.1489
2 -49365.4297 -52948.2459
3 -43460.1078 -49365.4297
4 -48294.6147 -43460.1078
5 -38682.2459 -48294.6147
6 -31278.0336 -38682.2459
7 -28067.9240 -31278.0336
8 -21399.2801 -28067.9240
9 -21169.2288 -21399.2801
10 -21382.9468 -21169.2288
11 -18723.7459 -21382.9468
12 -3480.9989 -18723.7459
13 -5600.8563 -3480.9989
14 -4474.9374 -5600.8563
15 -4148.7524 -4474.9374
16 1708.5011 -4148.7524
17 571.1437 1708.5011
18 -894.7124 571.1437
19 6483.4657 -894.7124
20 12713.2122 6483.4657
21 34877.9897 12713.2122
22 31616.3744 34877.9897
23 26781.6095 31616.3744
24 31235.3223 26781.6095
25 35785.2938 31235.3223
26 46840.9731 35785.2938
27 35539.4661 46840.9731
28 38678.8223 35539.4661
29 39738.2938 38678.8223
30 39171.4376 39738.2938
31 38802.5815 39171.4376
32 43867.2938 38802.5815
33 33825.5163 43867.2938
34 37692.7298 33825.5163
35 39992.7938 37692.7298
36 39954.6435 39992.7938
37 39980.6834 39954.6435
38 29670.7735 39980.6834
39 34798.9927 29670.7735
40 44826.2804 34798.9927
41 39132.8888 44826.2804
42 30810.1011 39132.8888
43 29871.2107 30810.1011
44 32368.8546 29871.2107
45 15509.8375 32368.8546
46 10196.1537 15509.8375
47 -1719.6797 10196.1537
48 10625.0331 -1719.6797
49 3602.0388 10625.0331
50 -6931.0423 3602.0388
51 -7310.6861 -6931.0423
52 -21430.1246 -7310.6861
53 -36992.5846 -21430.1246
54 -32566.4408 -36992.5846
55 -34529.2285 -32566.4408
56 -49072.5846 -34529.2285
57 -42879.7387 -49072.5846
58 -38186.5594 -42879.7387
59 -31874.5296 -38186.5594
60 -18830.8511 -31874.5296
61 -20818.9138 -18830.8511
62 -15740.3372 -20818.9138
63 -15418.9126 -15740.3372
64 -15488.8645 -15418.9126
65 -3767.4958 -15488.8645
66 -5242.3519 -3767.4958
67 -12560.1054 -5242.3519
68 -18477.4958 -12560.1054
69 -20164.3760 -18477.4958
70 -19935.7517 -20164.3760
71 -14456.4481 -19935.7517
> 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/7079b1258743587.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/82nrn1258743587.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/936ve1258743587.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/10u6zb1258743587.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/11y7pa1258743587.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/12n23a1258743587.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/13ubgg1258743587.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/14u13w1258743587.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/15nuc61258743587.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/16sk5m1258743587.tab")
+ }
>
> system("convert tmp/1177e1258743587.ps tmp/1177e1258743587.png")
> system("convert tmp/2p8tx1258743587.ps tmp/2p8tx1258743587.png")
> system("convert tmp/3fzi11258743587.ps tmp/3fzi11258743587.png")
> system("convert tmp/41ja81258743587.ps tmp/41ja81258743587.png")
> system("convert tmp/5n7fg1258743587.ps tmp/5n7fg1258743587.png")
> system("convert tmp/6h3o81258743587.ps tmp/6h3o81258743587.png")
> system("convert tmp/7079b1258743587.ps tmp/7079b1258743587.png")
> system("convert tmp/82nrn1258743587.ps tmp/82nrn1258743587.png")
> system("convert tmp/936ve1258743587.ps tmp/936ve1258743587.png")
> system("convert tmp/10u6zb1258743587.ps tmp/10u6zb1258743587.png")
>
>
> proc.time()
user system elapsed
2.439 1.554 2.962