R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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(492865,0,480961,0,461935,0,456608,0,441977,0,439148,0,488180,0,520564,0,501492,0,485025,0,464196,0,460170,0,467037,0,460070,0,447988,0,442867,0,436087,0,431328,0,484015,0,509673,0,512927,0,502831,0,470984,0,471067,0,476049,0,474605,0,470439,0,461251,0,454724,0,455626,0,516847,0,525192,0,522975,0,518585,0,509239,0,512238,0,519164,0,517009,0,509933,0,509127,0,500857,0,506971,0,569323,0,579714,0,577992,0,565464,0,547344,0,554788,0,562325,0,560854,0,555332,0,543599,0,536662,0,542722,1,593530,1,610763,1,612613,1,611324,1,594167,1,595454,1,590865,1,589379,1,584428,1,573100,1,567456,1,569028,1,620735,1,628884,1,628232,1,612117,1,595404,1,597141,1,593408,1,590072,1,579799,1,574205,1,572775,1,572942,1,619567,1,625809,1,619916,1,587625,1,565742,1,557274,1,560576,1,548854,0,531673,0,525919,0,511038,0,498662,0,555362,0,564591,0,541657,0,527070,0,509846,0,514258,0,516922,0,507561,0,492622,0,490243,0,469357,0,477580,0,528379,0,533590,0),dim=c(2,104),dimnames=list(c('W','D'),1:104))
> y <- array(NA,dim=c(2,104),dimnames=list(c('W','D'),1:104))
> 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
W D
1 492865 0
2 480961 0
3 461935 0
4 456608 0
5 441977 0
6 439148 0
7 488180 0
8 520564 0
9 501492 0
10 485025 0
11 464196 0
12 460170 0
13 467037 0
14 460070 0
15 447988 0
16 442867 0
17 436087 0
18 431328 0
19 484015 0
20 509673 0
21 512927 0
22 502831 0
23 470984 0
24 471067 0
25 476049 0
26 474605 0
27 470439 0
28 461251 0
29 454724 0
30 455626 0
31 516847 0
32 525192 0
33 522975 0
34 518585 0
35 509239 0
36 512238 0
37 519164 0
38 517009 0
39 509933 0
40 509127 0
41 500857 0
42 506971 0
43 569323 0
44 579714 0
45 577992 0
46 565464 0
47 547344 0
48 554788 0
49 562325 0
50 560854 0
51 555332 0
52 543599 0
53 536662 0
54 542722 1
55 593530 1
56 610763 1
57 612613 1
58 611324 1
59 594167 1
60 595454 1
61 590865 1
62 589379 1
63 584428 1
64 573100 1
65 567456 1
66 569028 1
67 620735 1
68 628884 1
69 628232 1
70 612117 1
71 595404 1
72 597141 1
73 593408 1
74 590072 1
75 579799 1
76 574205 1
77 572775 1
78 572942 1
79 619567 1
80 625809 1
81 619916 1
82 587625 1
83 565742 1
84 557274 1
85 560576 1
86 548854 0
87 531673 0
88 525919 0
89 511038 0
90 498662 0
91 555362 0
92 564591 0
93 541657 0
94 527070 0
95 509846 0
96 514258 0
97 516922 0
98 507561 0
99 492622 0
100 490243 0
101 469357 0
102 477580 0
103 528379 0
104 533590 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) D
504020 87763
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-72692 -24755 2668 21354 75694
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 504020 4010 125.69 <2e-16 ***
D 87763 7229 12.14 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 34030 on 102 degrees of freedom
Multiple R-squared: 0.591, Adjusted R-squared: 0.587
F-statistic: 147.4 on 1 and 102 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.3040113 0.6080225915 6.959887e-01
[2,] 0.2790810 0.5581620347 7.209190e-01
[3,] 0.2385420 0.4770840683 7.614580e-01
[4,] 0.4444021 0.8888042114 5.555979e-01
[5,] 0.4035300 0.8070599993 5.964700e-01
[6,] 0.3063742 0.6127484784 6.936258e-01
[7,] 0.2391375 0.4782750604 7.608625e-01
[8,] 0.1925686 0.3851371799 8.074314e-01
[9,] 0.1415591 0.2831182488 8.584409e-01
[10,] 0.1118107 0.2236214358 8.881893e-01
[11,] 0.1143360 0.2286719325 8.856640e-01
[12,] 0.1333100 0.2666200266 8.666900e-01
[13,] 0.1849599 0.3699198778 8.150401e-01
[14,] 0.2788307 0.5576613402 7.211693e-01
[15,] 0.2501599 0.5003198963 7.498401e-01
[16,] 0.3148269 0.6296537512 6.851731e-01
[17,] 0.3800541 0.7601081421 6.199459e-01
[18,] 0.3813471 0.7626941017 6.186529e-01
[19,] 0.3473500 0.6947000841 6.526500e-01
[20,] 0.3183442 0.6366883839 6.816558e-01
[21,] 0.2881961 0.5763922582 7.118039e-01
[22,] 0.2642345 0.5284690224 7.357655e-01
[23,] 0.2515239 0.5030478527 7.484761e-01
[24,] 0.2698318 0.5396635030 7.301682e-01
[25,] 0.3303265 0.6606529132 6.696735e-01
[26,] 0.4106264 0.8212527493 5.893736e-01
[27,] 0.5010979 0.9978042324 4.989021e-01
[28,] 0.6165421 0.7669158413 3.834579e-01
[29,] 0.6888642 0.6222716164 3.111358e-01
[30,] 0.7248912 0.5502175273 2.751088e-01
[31,] 0.7292886 0.5414228311 2.707114e-01
[32,] 0.7359934 0.5280131336 2.640066e-01
[33,] 0.7535487 0.4929026299 2.464513e-01
[34,] 0.7607571 0.4784857275 2.392429e-01
[35,] 0.7545833 0.4908333010 2.454167e-01
[36,] 0.7466015 0.5067970259 2.533985e-01
[37,] 0.7370268 0.5259464182 2.629732e-01
[38,] 0.7282284 0.5435431864 2.717716e-01
[39,] 0.8997013 0.2005973326 1.002987e-01
[40,] 0.9795336 0.0409328764 2.046644e-02
[41,] 0.9961141 0.0077718071 3.885904e-03
[42,] 0.9986108 0.0027783014 1.389151e-03
[43,] 0.9988899 0.0022201244 1.110062e-03
[44,] 0.9993085 0.0013829949 6.914975e-04
[45,] 0.9996959 0.0006082603 3.041302e-04
[46,] 0.9998595 0.0002810034 1.405017e-04
[47,] 0.9999162 0.0001676280 8.381399e-05
[48,] 0.9999176 0.0001647070 8.235348e-05
[49,] 0.9998987 0.0002025005 1.012503e-04
[50,] 0.9999465 0.0001070865 5.354325e-05
[51,] 0.9999198 0.0001603003 8.015016e-05
[52,] 0.9998967 0.0002066199 1.033099e-04
[53,] 0.9998633 0.0002733616 1.366808e-04
[54,] 0.9998096 0.0003807782 1.903891e-04
[55,] 0.9996680 0.0006639032 3.319516e-04
[56,] 0.9994348 0.0011303113 5.651556e-04
[57,] 0.9990488 0.0019023571 9.511785e-04
[58,] 0.9984329 0.0031342942 1.567147e-03
[59,] 0.9975257 0.0049485249 2.474262e-03
[60,] 0.9966835 0.0066329678 3.316484e-03
[61,] 0.9961269 0.0077461727 3.873086e-03
[62,] 0.9954027 0.0091946917 4.597346e-03
[63,] 0.9950505 0.0098989885 4.949494e-03
[64,] 0.9960466 0.0079068824 3.953441e-03
[65,] 0.9970111 0.0059777179 2.988859e-03
[66,] 0.9963573 0.0072854135 3.642707e-03
[67,] 0.9943003 0.0113993612 5.699681e-03
[68,] 0.9914074 0.0171851262 8.592563e-03
[69,] 0.9869691 0.0260617743 1.303089e-02
[70,] 0.9803809 0.0392381863 1.961909e-02
[71,] 0.9715280 0.0569439606 2.847198e-02
[72,] 0.9613481 0.0773037160 3.865186e-02
[73,] 0.9495228 0.1009543960 5.047720e-02
[74,] 0.9357079 0.1285841232 6.429206e-02
[75,] 0.9349641 0.1300717146 6.503586e-02
[76,] 0.9533173 0.0933653200 4.668266e-02
[77,] 0.9721123 0.0557753479 2.788767e-02
[78,] 0.9661934 0.0676131921 3.380660e-02
[79,] 0.9514145 0.0971709594 4.858548e-02
[80,] 0.9321859 0.1356281107 6.781406e-02
[81,] 0.9054880 0.1890240706 9.451204e-02
[82,] 0.9132835 0.1734330894 8.671654e-02
[83,] 0.8888005 0.2223990542 1.111995e-01
[84,] 0.8505529 0.2988942833 1.494471e-01
[85,] 0.7940956 0.4118087506 2.059044e-01
[86,] 0.7425826 0.5148347428 2.574174e-01
[87,] 0.7946693 0.4106614390 2.053307e-01
[88,] 0.9152945 0.1694109744 8.470549e-02
[89,] 0.9322445 0.1355109708 6.775549e-02
[90,] 0.9157800 0.1684400394 8.422002e-02
[91,] 0.8596111 0.2807778389 1.403889e-01
[92,] 0.7872348 0.4255304019 2.127652e-01
[93,] 0.7037924 0.5924152961 2.962076e-01
[94,] 0.5688623 0.8622754467 4.311377e-01
[95,] 0.4056528 0.8113056010 5.943472e-01
> postscript(file="/var/www/html/freestat/rcomp/tmp/1sizo1227790159.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/freestat/rcomp/tmp/2d83y1227790159.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/freestat/rcomp/tmp/3pj4s1227790159.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/freestat/rcomp/tmp/4elnv1227790159.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/freestat/rcomp/tmp/5bphk1227790159.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 = 104
Frequency = 1
1 2 3 4 5 6 7
-11154.542 -23058.542 -42084.542 -47411.542 -62042.542 -64871.542 -15839.542
8 9 10 11 12 13 14
16544.458 -2527.542 -18994.542 -39823.542 -43849.542 -36982.542 -43949.542
15 16 17 18 19 20 21
-56031.542 -61152.542 -67932.542 -72691.542 -20004.542 5653.458 8907.458
22 23 24 25 26 27 28
-1188.542 -33035.542 -32952.542 -27970.542 -29414.542 -33580.542 -42768.542
29 30 31 32 33 34 35
-49295.542 -48393.542 12827.458 21172.458 18955.458 14565.458 5219.458
36 37 38 39 40 41 42
8218.458 15144.458 12989.458 5913.458 5107.458 -3162.542 2951.458
43 44 45 46 47 48 49
65303.458 75694.458 73972.458 61444.458 43324.458 50768.458 58305.458
50 51 52 53 54 55 56
56834.458 51312.458 39579.458 32642.458 -49060.875 1747.125 18980.125
57 58 59 60 61 62 63
20830.125 19541.125 2384.125 3671.125 -917.875 -2403.875 -7354.875
64 65 66 67 68 69 70
-18682.875 -24326.875 -22754.875 28952.125 37101.125 36449.125 20334.125
71 72 73 74 75 76 77
3621.125 5358.125 1625.125 -1710.875 -11983.875 -17577.875 -19007.875
78 79 80 81 82 83 84
-18840.875 27784.125 34026.125 28133.125 -4157.875 -26040.875 -34508.875
85 86 87 88 89 90 91
-31206.875 44834.458 27653.458 21899.458 7018.458 -5357.542 51342.458
92 93 94 95 96 97 98
60571.458 37637.458 23050.458 5826.458 10238.458 12902.458 3541.458
99 100 101 102 103 104
-11397.542 -13776.542 -34662.542 -26439.542 24359.458 29570.458
> postscript(file="/var/www/html/freestat/rcomp/tmp/6hj1z1227790159.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 = 104
Frequency = 1
lag(myerror, k = 1) myerror
0 -11154.542 NA
1 -23058.542 -11154.542
2 -42084.542 -23058.542
3 -47411.542 -42084.542
4 -62042.542 -47411.542
5 -64871.542 -62042.542
6 -15839.542 -64871.542
7 16544.458 -15839.542
8 -2527.542 16544.458
9 -18994.542 -2527.542
10 -39823.542 -18994.542
11 -43849.542 -39823.542
12 -36982.542 -43849.542
13 -43949.542 -36982.542
14 -56031.542 -43949.542
15 -61152.542 -56031.542
16 -67932.542 -61152.542
17 -72691.542 -67932.542
18 -20004.542 -72691.542
19 5653.458 -20004.542
20 8907.458 5653.458
21 -1188.542 8907.458
22 -33035.542 -1188.542
23 -32952.542 -33035.542
24 -27970.542 -32952.542
25 -29414.542 -27970.542
26 -33580.542 -29414.542
27 -42768.542 -33580.542
28 -49295.542 -42768.542
29 -48393.542 -49295.542
30 12827.458 -48393.542
31 21172.458 12827.458
32 18955.458 21172.458
33 14565.458 18955.458
34 5219.458 14565.458
35 8218.458 5219.458
36 15144.458 8218.458
37 12989.458 15144.458
38 5913.458 12989.458
39 5107.458 5913.458
40 -3162.542 5107.458
41 2951.458 -3162.542
42 65303.458 2951.458
43 75694.458 65303.458
44 73972.458 75694.458
45 61444.458 73972.458
46 43324.458 61444.458
47 50768.458 43324.458
48 58305.458 50768.458
49 56834.458 58305.458
50 51312.458 56834.458
51 39579.458 51312.458
52 32642.458 39579.458
53 -49060.875 32642.458
54 1747.125 -49060.875
55 18980.125 1747.125
56 20830.125 18980.125
57 19541.125 20830.125
58 2384.125 19541.125
59 3671.125 2384.125
60 -917.875 3671.125
61 -2403.875 -917.875
62 -7354.875 -2403.875
63 -18682.875 -7354.875
64 -24326.875 -18682.875
65 -22754.875 -24326.875
66 28952.125 -22754.875
67 37101.125 28952.125
68 36449.125 37101.125
69 20334.125 36449.125
70 3621.125 20334.125
71 5358.125 3621.125
72 1625.125 5358.125
73 -1710.875 1625.125
74 -11983.875 -1710.875
75 -17577.875 -11983.875
76 -19007.875 -17577.875
77 -18840.875 -19007.875
78 27784.125 -18840.875
79 34026.125 27784.125
80 28133.125 34026.125
81 -4157.875 28133.125
82 -26040.875 -4157.875
83 -34508.875 -26040.875
84 -31206.875 -34508.875
85 44834.458 -31206.875
86 27653.458 44834.458
87 21899.458 27653.458
88 7018.458 21899.458
89 -5357.542 7018.458
90 51342.458 -5357.542
91 60571.458 51342.458
92 37637.458 60571.458
93 23050.458 37637.458
94 5826.458 23050.458
95 10238.458 5826.458
96 12902.458 10238.458
97 3541.458 12902.458
98 -11397.542 3541.458
99 -13776.542 -11397.542
100 -34662.542 -13776.542
101 -26439.542 -34662.542
102 24359.458 -26439.542
103 29570.458 24359.458
104 NA 29570.458
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -23058.542 -11154.542
[2,] -42084.542 -23058.542
[3,] -47411.542 -42084.542
[4,] -62042.542 -47411.542
[5,] -64871.542 -62042.542
[6,] -15839.542 -64871.542
[7,] 16544.458 -15839.542
[8,] -2527.542 16544.458
[9,] -18994.542 -2527.542
[10,] -39823.542 -18994.542
[11,] -43849.542 -39823.542
[12,] -36982.542 -43849.542
[13,] -43949.542 -36982.542
[14,] -56031.542 -43949.542
[15,] -61152.542 -56031.542
[16,] -67932.542 -61152.542
[17,] -72691.542 -67932.542
[18,] -20004.542 -72691.542
[19,] 5653.458 -20004.542
[20,] 8907.458 5653.458
[21,] -1188.542 8907.458
[22,] -33035.542 -1188.542
[23,] -32952.542 -33035.542
[24,] -27970.542 -32952.542
[25,] -29414.542 -27970.542
[26,] -33580.542 -29414.542
[27,] -42768.542 -33580.542
[28,] -49295.542 -42768.542
[29,] -48393.542 -49295.542
[30,] 12827.458 -48393.542
[31,] 21172.458 12827.458
[32,] 18955.458 21172.458
[33,] 14565.458 18955.458
[34,] 5219.458 14565.458
[35,] 8218.458 5219.458
[36,] 15144.458 8218.458
[37,] 12989.458 15144.458
[38,] 5913.458 12989.458
[39,] 5107.458 5913.458
[40,] -3162.542 5107.458
[41,] 2951.458 -3162.542
[42,] 65303.458 2951.458
[43,] 75694.458 65303.458
[44,] 73972.458 75694.458
[45,] 61444.458 73972.458
[46,] 43324.458 61444.458
[47,] 50768.458 43324.458
[48,] 58305.458 50768.458
[49,] 56834.458 58305.458
[50,] 51312.458 56834.458
[51,] 39579.458 51312.458
[52,] 32642.458 39579.458
[53,] -49060.875 32642.458
[54,] 1747.125 -49060.875
[55,] 18980.125 1747.125
[56,] 20830.125 18980.125
[57,] 19541.125 20830.125
[58,] 2384.125 19541.125
[59,] 3671.125 2384.125
[60,] -917.875 3671.125
[61,] -2403.875 -917.875
[62,] -7354.875 -2403.875
[63,] -18682.875 -7354.875
[64,] -24326.875 -18682.875
[65,] -22754.875 -24326.875
[66,] 28952.125 -22754.875
[67,] 37101.125 28952.125
[68,] 36449.125 37101.125
[69,] 20334.125 36449.125
[70,] 3621.125 20334.125
[71,] 5358.125 3621.125
[72,] 1625.125 5358.125
[73,] -1710.875 1625.125
[74,] -11983.875 -1710.875
[75,] -17577.875 -11983.875
[76,] -19007.875 -17577.875
[77,] -18840.875 -19007.875
[78,] 27784.125 -18840.875
[79,] 34026.125 27784.125
[80,] 28133.125 34026.125
[81,] -4157.875 28133.125
[82,] -26040.875 -4157.875
[83,] -34508.875 -26040.875
[84,] -31206.875 -34508.875
[85,] 44834.458 -31206.875
[86,] 27653.458 44834.458
[87,] 21899.458 27653.458
[88,] 7018.458 21899.458
[89,] -5357.542 7018.458
[90,] 51342.458 -5357.542
[91,] 60571.458 51342.458
[92,] 37637.458 60571.458
[93,] 23050.458 37637.458
[94,] 5826.458 23050.458
[95,] 10238.458 5826.458
[96,] 12902.458 10238.458
[97,] 3541.458 12902.458
[98,] -11397.542 3541.458
[99,] -13776.542 -11397.542
[100,] -34662.542 -13776.542
[101,] -26439.542 -34662.542
[102,] 24359.458 -26439.542
[103,] 29570.458 24359.458
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -23058.542 -11154.542
2 -42084.542 -23058.542
3 -47411.542 -42084.542
4 -62042.542 -47411.542
5 -64871.542 -62042.542
6 -15839.542 -64871.542
7 16544.458 -15839.542
8 -2527.542 16544.458
9 -18994.542 -2527.542
10 -39823.542 -18994.542
11 -43849.542 -39823.542
12 -36982.542 -43849.542
13 -43949.542 -36982.542
14 -56031.542 -43949.542
15 -61152.542 -56031.542
16 -67932.542 -61152.542
17 -72691.542 -67932.542
18 -20004.542 -72691.542
19 5653.458 -20004.542
20 8907.458 5653.458
21 -1188.542 8907.458
22 -33035.542 -1188.542
23 -32952.542 -33035.542
24 -27970.542 -32952.542
25 -29414.542 -27970.542
26 -33580.542 -29414.542
27 -42768.542 -33580.542
28 -49295.542 -42768.542
29 -48393.542 -49295.542
30 12827.458 -48393.542
31 21172.458 12827.458
32 18955.458 21172.458
33 14565.458 18955.458
34 5219.458 14565.458
35 8218.458 5219.458
36 15144.458 8218.458
37 12989.458 15144.458
38 5913.458 12989.458
39 5107.458 5913.458
40 -3162.542 5107.458
41 2951.458 -3162.542
42 65303.458 2951.458
43 75694.458 65303.458
44 73972.458 75694.458
45 61444.458 73972.458
46 43324.458 61444.458
47 50768.458 43324.458
48 58305.458 50768.458
49 56834.458 58305.458
50 51312.458 56834.458
51 39579.458 51312.458
52 32642.458 39579.458
53 -49060.875 32642.458
54 1747.125 -49060.875
55 18980.125 1747.125
56 20830.125 18980.125
57 19541.125 20830.125
58 2384.125 19541.125
59 3671.125 2384.125
60 -917.875 3671.125
61 -2403.875 -917.875
62 -7354.875 -2403.875
63 -18682.875 -7354.875
64 -24326.875 -18682.875
65 -22754.875 -24326.875
66 28952.125 -22754.875
67 37101.125 28952.125
68 36449.125 37101.125
69 20334.125 36449.125
70 3621.125 20334.125
71 5358.125 3621.125
72 1625.125 5358.125
73 -1710.875 1625.125
74 -11983.875 -1710.875
75 -17577.875 -11983.875
76 -19007.875 -17577.875
77 -18840.875 -19007.875
78 27784.125 -18840.875
79 34026.125 27784.125
80 28133.125 34026.125
81 -4157.875 28133.125
82 -26040.875 -4157.875
83 -34508.875 -26040.875
84 -31206.875 -34508.875
85 44834.458 -31206.875
86 27653.458 44834.458
87 21899.458 27653.458
88 7018.458 21899.458
89 -5357.542 7018.458
90 51342.458 -5357.542
91 60571.458 51342.458
92 37637.458 60571.458
93 23050.458 37637.458
94 5826.458 23050.458
95 10238.458 5826.458
96 12902.458 10238.458
97 3541.458 12902.458
98 -11397.542 3541.458
99 -13776.542 -11397.542
100 -34662.542 -13776.542
101 -26439.542 -34662.542
102 24359.458 -26439.542
103 29570.458 24359.458
> 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/freestat/rcomp/tmp/7f3c61227790159.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/freestat/rcomp/tmp/8jr6z1227790160.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/freestat/rcomp/tmp/916ez1227790160.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/freestat/rcomp/tmp/10gpmf1227790160.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11lwnu1227790160.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/freestat/rcomp/tmp/12blq61227790160.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/freestat/rcomp/tmp/13pnf81227790160.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/freestat/rcomp/tmp/1428291227790160.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/freestat/rcomp/tmp/153f0s1227790160.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/freestat/rcomp/tmp/163iui1227790160.tab")
+ }
>
> system("convert tmp/1sizo1227790159.ps tmp/1sizo1227790159.png")
> system("convert tmp/2d83y1227790159.ps tmp/2d83y1227790159.png")
> system("convert tmp/3pj4s1227790159.ps tmp/3pj4s1227790159.png")
> system("convert tmp/4elnv1227790159.ps tmp/4elnv1227790159.png")
> system("convert tmp/5bphk1227790159.ps tmp/5bphk1227790159.png")
> system("convert tmp/6hj1z1227790159.ps tmp/6hj1z1227790159.png")
> system("convert tmp/7f3c61227790159.ps tmp/7f3c61227790159.png")
> system("convert tmp/8jr6z1227790160.ps tmp/8jr6z1227790160.png")
> system("convert tmp/916ez1227790160.ps tmp/916ez1227790160.png")
> system("convert tmp/10gpmf1227790160.ps tmp/10gpmf1227790160.png")
>
>
> proc.time()
user system elapsed
4.122 2.476 4.513