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 = 'No 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
1 519164 0.9 1 0 0 0 0 0 0 0 0 0 0
2 517009 1.3 0 1 0 0 0 0 0 0 0 0 0
3 509933 1.4 0 0 1 0 0 0 0 0 0 0 0
4 509127 1.5 0 0 0 1 0 0 0 0 0 0 0
5 500857 1.1 0 0 0 0 1 0 0 0 0 0 0
6 506971 1.6 0 0 0 0 0 1 0 0 0 0 0
7 569323 1.5 0 0 0 0 0 0 1 0 0 0 0
8 579714 1.6 0 0 0 0 0 0 0 1 0 0 0
9 577992 1.7 0 0 0 0 0 0 0 0 1 0 0
10 565464 1.6 0 0 0 0 0 0 0 0 0 1 0
11 547344 1.7 0 0 0 0 0 0 0 0 0 0 1
12 554788 1.6 0 0 0 0 0 0 0 0 0 0 0
13 562325 1.6 1 0 0 0 0 0 0 0 0 0 0
14 560854 1.3 0 1 0 0 0 0 0 0 0 0 0
15 555332 1.1 0 0 1 0 0 0 0 0 0 0 0
16 543599 1.6 0 0 0 1 0 0 0 0 0 0 0
17 536662 1.9 0 0 0 0 1 0 0 0 0 0 0
18 542722 1.6 0 0 0 0 0 1 0 0 0 0 0
19 593530 1.7 0 0 0 0 0 0 1 0 0 0 0
20 610763 1.6 0 0 0 0 0 0 0 1 0 0 0
21 612613 1.4 0 0 0 0 0 0 0 0 1 0 0
22 611324 2.1 0 0 0 0 0 0 0 0 0 1 0
23 594167 1.9 0 0 0 0 0 0 0 0 0 0 1
24 595454 1.7 0 0 0 0 0 0 0 0 0 0 0
25 590865 1.8 1 0 0 0 0 0 0 0 0 0 0
26 589379 2.0 0 1 0 0 0 0 0 0 0 0 0
27 584428 2.5 0 0 1 0 0 0 0 0 0 0 0
28 573100 2.1 0 0 0 1 0 0 0 0 0 0 0
29 567456 2.1 0 0 0 0 1 0 0 0 0 0 0
30 569028 2.3 0 0 0 0 0 1 0 0 0 0 0
31 620735 2.4 0 0 0 0 0 0 1 0 0 0 0
32 628884 2.4 0 0 0 0 0 0 0 1 0 0 0
33 628232 2.3 0 0 0 0 0 0 0 0 1 0 0
34 612117 1.7 0 0 0 0 0 0 0 0 0 1 0
35 595404 2.0 0 0 0 0 0 0 0 0 0 0 1
36 597141 2.3 0 0 0 0 0 0 0 0 0 0 0
37 593408 2.0 1 0 0 0 0 0 0 0 0 0 0
38 590072 2.0 0 1 0 0 0 0 0 0 0 0 0
39 579799 1.3 0 0 1 0 0 0 0 0 0 0 0
40 574205 1.7 0 0 0 1 0 0 0 0 0 0 0
41 572775 1.9 0 0 0 0 1 0 0 0 0 0 0
42 572942 1.7 0 0 0 0 0 1 0 0 0 0 0
43 619567 1.6 0 0 0 0 0 0 1 0 0 0 0
44 625809 1.7 0 0 0 0 0 0 0 1 0 0 0
45 619916 1.8 0 0 0 0 0 0 0 0 1 0 0
46 587625 1.9 0 0 0 0 0 0 0 0 0 1 0
47 565742 1.9 0 0 0 0 0 0 0 0 0 0 1
48 557274 1.9 0 0 0 0 0 0 0 0 0 0 0
49 560576 2.0 1 0 0 0 0 0 0 0 0 0 0
50 548854 2.1 0 1 0 0 0 0 0 0 0 0 0
51 531673 1.9 0 0 1 0 0 0 0 0 0 0 0
52 525919 1.9 0 0 0 1 0 0 0 0 0 0 0
53 511038 1.3 0 0 0 0 1 0 0 0 0 0 0
54 498662 1.3 0 0 0 0 0 1 0 0 0 0 0
55 555362 1.4 0 0 0 0 0 0 1 0 0 0 0
56 564591 1.2 0 0 0 0 0 0 0 1 0 0 0
57 541657 1.3 0 0 0 0 0 0 0 0 1 0 0
58 527070 1.8 0 0 0 0 0 0 0 0 0 1 0
59 509846 2.2 0 0 0 0 0 0 0 0 0 0 1
60 514258 2.6 0 0 0 0 0 0 0 0 0 0 0
61 516922 2.8 1 0 0 0 0 0 0 0 0 0 0
62 507561 3.1 0 1 0 0 0 0 0 0 0 0 0
63 492622 3.9 0 0 1 0 0 0 0 0 0 0 0
64 490243 3.7 0 0 0 1 0 0 0 0 0 0 0
65 469357 4.6 0 0 0 0 1 0 0 0 0 0 0
66 477580 5.1 0 0 0 0 0 1 0 0 0 0 0
67 528379 5.2 0 0 0 0 0 0 1 0 0 0 0
68 533590 4.9 0 0 0 0 0 0 0 1 0 0 0
69 517945 5.1 0 0 0 0 0 0 0 0 1 0 0
70 506174 4.8 0 0 0 0 0 0 0 0 0 1 0
71 501866 3.9 0 0 0 0 0 0 0 0 0 0 1
72 516141 3.5 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) GI M1 M2 M3 M4
594928 -17244 -5817 -8728 -17856 -22972
M5 M6 M7 M8 M9 M10
-31497 -27858 25881 34141 27216 13315
M11
-3448
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-58071 -21915 -2156 31870 50465
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 594928 16315 36.465 < 2e-16 ***
GI -17244 3888 -4.435 4.08e-05 ***
M1 -5818 19484 -0.299 0.7663
M2 -8728 19452 -0.449 0.6553
M3 -17856 19441 -0.918 0.3621
M4 -22972 19430 -1.182 0.2418
M5 -31497 19422 -1.622 0.1102
M6 -27858 19417 -1.435 0.1566
M7 25882 19417 1.333 0.1877
M8 34141 19417 1.758 0.0839 .
M9 27216 19417 1.402 0.1662
M10 13315 19418 0.686 0.4956
M11 -3448 19417 -0.178 0.8597
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 33630 on 59 degrees of freedom
Multiple R-squared: 0.4388, Adjusted R-squared: 0.3246
F-statistic: 3.844 on 12 and 59 DF, p-value: 0.0002426
> 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.572496108 0.8550077849 4.275039e-01
[2,] 0.401270357 0.8025407142 5.987296e-01
[3,] 0.339616141 0.6792322830 6.603839e-01
[4,] 0.233073132 0.4661462647 7.669269e-01
[5,] 0.182826862 0.3656537245 8.171731e-01
[6,] 0.183636721 0.3672734416 8.163633e-01
[7,] 0.148430734 0.2968614689 8.515693e-01
[8,] 0.137543387 0.2750867748 8.624566e-01
[9,] 0.120425357 0.2408507131 8.795746e-01
[10,] 0.089255466 0.1785109328 9.107445e-01
[11,] 0.059954723 0.1199094452 9.400453e-01
[12,] 0.048135030 0.0962700603 9.518650e-01
[13,] 0.034858391 0.0697167825 9.651416e-01
[14,] 0.024990127 0.0499802539 9.750099e-01
[15,] 0.017097334 0.0341946673 9.829027e-01
[16,] 0.011998524 0.0239970473 9.880015e-01
[17,] 0.009210802 0.0184216044 9.907892e-01
[18,] 0.008092315 0.0161846308 9.919077e-01
[19,] 0.009721925 0.0194438505 9.902781e-01
[20,] 0.009545157 0.0190903145 9.904548e-01
[21,] 0.010712019 0.0214240374 9.892880e-01
[22,] 0.008855597 0.0177111949 9.911444e-01
[23,] 0.008258026 0.0165160529 9.917420e-01
[24,] 0.018144916 0.0362898324 9.818551e-01
[25,] 0.025001609 0.0500032175 9.749984e-01
[26,] 0.045419858 0.0908397159 9.545801e-01
[27,] 0.097041000 0.1940820004 9.029590e-01
[28,] 0.159283683 0.3185673668 8.407163e-01
[29,] 0.258854817 0.5177096343 7.411452e-01
[30,] 0.559551158 0.8808976846 4.404488e-01
[31,] 0.778806288 0.4423874245 2.211937e-01
[32,] 0.923957390 0.1520852207 7.604261e-02
[33,] 0.957762323 0.0844753544 4.223768e-02
[34,] 0.984318469 0.0313630613 1.568153e-02
[35,] 0.996211872 0.0075762557 3.788128e-03
[36,] 0.998248770 0.0035024605 1.751230e-03
[37,] 0.999485148 0.0010297049 5.148525e-04
[38,] 0.999943049 0.0001139015 5.695076e-05
[39,] 0.999754196 0.0004916081 2.458040e-04
[40,] 0.998388509 0.0032229826 1.611491e-03
[41,] 0.995462379 0.0090752424 4.537621e-03
> postscript(file="/var/www/html/rcomp/tmp/1bwta1258743427.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/2f5r11258743427.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/3kveg1258743427.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/4em4y1258743427.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/5h01p1258743427.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
-54427.3406 -46774.8443 -42998.3351 -36963.8846 -43606.1923 -32508.8443
7 8 9 10 11 12
-25621.1465 -21765.3938 -14838.4927 -15189.5201 -14822.1593 -12550.3443
13 14 15 16 17 18
804.1209 -2929.8443 -2772.3901 -767.5330 5993.6209 3242.1557
19 20 21 22 23 24
2034.5568 9283.6062 14609.4524 39292.2381 35449.5440 29840.0073
25 26 27 28 29 30
32792.8242 37665.6172 50464.5330 37355.2253 40236.3242 41618.6172
31 32 33 34 35 36
41310.0183 41199.4194 45747.6172 33187.8315 38410.8956 41873.1172
37 38 39 40 41 42
38784.5275 38358.6172 25143.3132 31562.8187 42106.6209 35186.5073
43 44 45 46 47 48
26347.2052 26053.9579 28809.8590 12144.5348 7024.5440 -4891.2894
49 50 51 52 53 54
5952.5275 -1135.0311 -12636.5769 -13274.4780 -29976.4890 -45990.8992
55 56 57 58 59 60
-41306.4981 -43785.8003 -58070.8992 -50134.8168 -43698.4011 -35836.8278
61 62 63 64 65 66
-23906.6594 -25184.5147 -17200.5440 -17912.1484 -14753.8847 -1547.5367
67 68 69 70 71 72
-2764.1356 -10985.7895 -16257.5367 -19300.2675 -22364.4231 -18434.6630
> postscript(file="/var/www/html/rcomp/tmp/6y2nc1258743427.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 -54427.3406 NA
1 -46774.8443 -54427.3406
2 -42998.3351 -46774.8443
3 -36963.8846 -42998.3351
4 -43606.1923 -36963.8846
5 -32508.8443 -43606.1923
6 -25621.1465 -32508.8443
7 -21765.3938 -25621.1465
8 -14838.4927 -21765.3938
9 -15189.5201 -14838.4927
10 -14822.1593 -15189.5201
11 -12550.3443 -14822.1593
12 804.1209 -12550.3443
13 -2929.8443 804.1209
14 -2772.3901 -2929.8443
15 -767.5330 -2772.3901
16 5993.6209 -767.5330
17 3242.1557 5993.6209
18 2034.5568 3242.1557
19 9283.6062 2034.5568
20 14609.4524 9283.6062
21 39292.2381 14609.4524
22 35449.5440 39292.2381
23 29840.0073 35449.5440
24 32792.8242 29840.0073
25 37665.6172 32792.8242
26 50464.5330 37665.6172
27 37355.2253 50464.5330
28 40236.3242 37355.2253
29 41618.6172 40236.3242
30 41310.0183 41618.6172
31 41199.4194 41310.0183
32 45747.6172 41199.4194
33 33187.8315 45747.6172
34 38410.8956 33187.8315
35 41873.1172 38410.8956
36 38784.5275 41873.1172
37 38358.6172 38784.5275
38 25143.3132 38358.6172
39 31562.8187 25143.3132
40 42106.6209 31562.8187
41 35186.5073 42106.6209
42 26347.2052 35186.5073
43 26053.9579 26347.2052
44 28809.8590 26053.9579
45 12144.5348 28809.8590
46 7024.5440 12144.5348
47 -4891.2894 7024.5440
48 5952.5275 -4891.2894
49 -1135.0311 5952.5275
50 -12636.5769 -1135.0311
51 -13274.4780 -12636.5769
52 -29976.4890 -13274.4780
53 -45990.8992 -29976.4890
54 -41306.4981 -45990.8992
55 -43785.8003 -41306.4981
56 -58070.8992 -43785.8003
57 -50134.8168 -58070.8992
58 -43698.4011 -50134.8168
59 -35836.8278 -43698.4011
60 -23906.6594 -35836.8278
61 -25184.5147 -23906.6594
62 -17200.5440 -25184.5147
63 -17912.1484 -17200.5440
64 -14753.8847 -17912.1484
65 -1547.5367 -14753.8847
66 -2764.1356 -1547.5367
67 -10985.7895 -2764.1356
68 -16257.5367 -10985.7895
69 -19300.2675 -16257.5367
70 -22364.4231 -19300.2675
71 -18434.6630 -22364.4231
72 NA -18434.6630
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -46774.8443 -54427.3406
[2,] -42998.3351 -46774.8443
[3,] -36963.8846 -42998.3351
[4,] -43606.1923 -36963.8846
[5,] -32508.8443 -43606.1923
[6,] -25621.1465 -32508.8443
[7,] -21765.3938 -25621.1465
[8,] -14838.4927 -21765.3938
[9,] -15189.5201 -14838.4927
[10,] -14822.1593 -15189.5201
[11,] -12550.3443 -14822.1593
[12,] 804.1209 -12550.3443
[13,] -2929.8443 804.1209
[14,] -2772.3901 -2929.8443
[15,] -767.5330 -2772.3901
[16,] 5993.6209 -767.5330
[17,] 3242.1557 5993.6209
[18,] 2034.5568 3242.1557
[19,] 9283.6062 2034.5568
[20,] 14609.4524 9283.6062
[21,] 39292.2381 14609.4524
[22,] 35449.5440 39292.2381
[23,] 29840.0073 35449.5440
[24,] 32792.8242 29840.0073
[25,] 37665.6172 32792.8242
[26,] 50464.5330 37665.6172
[27,] 37355.2253 50464.5330
[28,] 40236.3242 37355.2253
[29,] 41618.6172 40236.3242
[30,] 41310.0183 41618.6172
[31,] 41199.4194 41310.0183
[32,] 45747.6172 41199.4194
[33,] 33187.8315 45747.6172
[34,] 38410.8956 33187.8315
[35,] 41873.1172 38410.8956
[36,] 38784.5275 41873.1172
[37,] 38358.6172 38784.5275
[38,] 25143.3132 38358.6172
[39,] 31562.8187 25143.3132
[40,] 42106.6209 31562.8187
[41,] 35186.5073 42106.6209
[42,] 26347.2052 35186.5073
[43,] 26053.9579 26347.2052
[44,] 28809.8590 26053.9579
[45,] 12144.5348 28809.8590
[46,] 7024.5440 12144.5348
[47,] -4891.2894 7024.5440
[48,] 5952.5275 -4891.2894
[49,] -1135.0311 5952.5275
[50,] -12636.5769 -1135.0311
[51,] -13274.4780 -12636.5769
[52,] -29976.4890 -13274.4780
[53,] -45990.8992 -29976.4890
[54,] -41306.4981 -45990.8992
[55,] -43785.8003 -41306.4981
[56,] -58070.8992 -43785.8003
[57,] -50134.8168 -58070.8992
[58,] -43698.4011 -50134.8168
[59,] -35836.8278 -43698.4011
[60,] -23906.6594 -35836.8278
[61,] -25184.5147 -23906.6594
[62,] -17200.5440 -25184.5147
[63,] -17912.1484 -17200.5440
[64,] -14753.8847 -17912.1484
[65,] -1547.5367 -14753.8847
[66,] -2764.1356 -1547.5367
[67,] -10985.7895 -2764.1356
[68,] -16257.5367 -10985.7895
[69,] -19300.2675 -16257.5367
[70,] -22364.4231 -19300.2675
[71,] -18434.6630 -22364.4231
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -46774.8443 -54427.3406
2 -42998.3351 -46774.8443
3 -36963.8846 -42998.3351
4 -43606.1923 -36963.8846
5 -32508.8443 -43606.1923
6 -25621.1465 -32508.8443
7 -21765.3938 -25621.1465
8 -14838.4927 -21765.3938
9 -15189.5201 -14838.4927
10 -14822.1593 -15189.5201
11 -12550.3443 -14822.1593
12 804.1209 -12550.3443
13 -2929.8443 804.1209
14 -2772.3901 -2929.8443
15 -767.5330 -2772.3901
16 5993.6209 -767.5330
17 3242.1557 5993.6209
18 2034.5568 3242.1557
19 9283.6062 2034.5568
20 14609.4524 9283.6062
21 39292.2381 14609.4524
22 35449.5440 39292.2381
23 29840.0073 35449.5440
24 32792.8242 29840.0073
25 37665.6172 32792.8242
26 50464.5330 37665.6172
27 37355.2253 50464.5330
28 40236.3242 37355.2253
29 41618.6172 40236.3242
30 41310.0183 41618.6172
31 41199.4194 41310.0183
32 45747.6172 41199.4194
33 33187.8315 45747.6172
34 38410.8956 33187.8315
35 41873.1172 38410.8956
36 38784.5275 41873.1172
37 38358.6172 38784.5275
38 25143.3132 38358.6172
39 31562.8187 25143.3132
40 42106.6209 31562.8187
41 35186.5073 42106.6209
42 26347.2052 35186.5073
43 26053.9579 26347.2052
44 28809.8590 26053.9579
45 12144.5348 28809.8590
46 7024.5440 12144.5348
47 -4891.2894 7024.5440
48 5952.5275 -4891.2894
49 -1135.0311 5952.5275
50 -12636.5769 -1135.0311
51 -13274.4780 -12636.5769
52 -29976.4890 -13274.4780
53 -45990.8992 -29976.4890
54 -41306.4981 -45990.8992
55 -43785.8003 -41306.4981
56 -58070.8992 -43785.8003
57 -50134.8168 -58070.8992
58 -43698.4011 -50134.8168
59 -35836.8278 -43698.4011
60 -23906.6594 -35836.8278
61 -25184.5147 -23906.6594
62 -17200.5440 -25184.5147
63 -17912.1484 -17200.5440
64 -14753.8847 -17912.1484
65 -1547.5367 -14753.8847
66 -2764.1356 -1547.5367
67 -10985.7895 -2764.1356
68 -16257.5367 -10985.7895
69 -19300.2675 -16257.5367
70 -22364.4231 -19300.2675
71 -18434.6630 -22364.4231
> 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/7shbq1258743427.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/8go2c1258743427.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/9fz2t1258743427.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/10v9841258743427.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/11qo4f1258743427.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/12dyi51258743427.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/1331b11258743427.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/14pz4s1258743427.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/1520gm1258743427.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/16ml9j1258743427.tab")
+ }
>
> system("convert tmp/1bwta1258743427.ps tmp/1bwta1258743427.png")
> system("convert tmp/2f5r11258743427.ps tmp/2f5r11258743427.png")
> system("convert tmp/3kveg1258743427.ps tmp/3kveg1258743427.png")
> system("convert tmp/4em4y1258743427.ps tmp/4em4y1258743427.png")
> system("convert tmp/5h01p1258743427.ps tmp/5h01p1258743427.png")
> system("convert tmp/6y2nc1258743427.ps tmp/6y2nc1258743427.png")
> system("convert tmp/7shbq1258743427.ps tmp/7shbq1258743427.png")
> system("convert tmp/8go2c1258743427.ps tmp/8go2c1258743427.png")
> system("convert tmp/9fz2t1258743427.ps tmp/9fz2t1258743427.png")
> system("convert tmp/10v9841258743427.ps tmp/10v9841258743427.png")
>
>
> proc.time()
user system elapsed
2.555 1.572 2.938