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,0,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 = '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 t
1 492865 0 1
2 480961 0 2
3 461935 0 3
4 456608 0 4
5 441977 0 5
6 439148 0 6
7 488180 0 7
8 520564 0 8
9 501492 0 9
10 485025 0 10
11 464196 0 11
12 460170 0 12
13 467037 0 13
14 460070 0 14
15 447988 0 15
16 442867 0 16
17 436087 0 17
18 431328 0 18
19 484015 0 19
20 509673 0 20
21 512927 0 21
22 502831 0 22
23 470984 0 23
24 471067 0 24
25 476049 0 25
26 474605 0 26
27 470439 0 27
28 461251 0 28
29 454724 0 29
30 455626 0 30
31 516847 0 31
32 525192 0 32
33 522975 0 33
34 518585 0 34
35 509239 0 35
36 512238 0 36
37 519164 0 37
38 517009 0 38
39 509933 0 39
40 509127 0 40
41 500857 0 41
42 506971 0 42
43 569323 0 43
44 579714 0 44
45 577992 0 45
46 565464 0 46
47 547344 0 47
48 554788 0 48
49 562325 0 49
50 560854 0 50
51 555332 0 51
52 543599 0 52
53 536662 0 53
54 542722 0 54
55 593530 1 55
56 610763 1 56
57 612613 1 57
58 611324 1 58
59 594167 1 59
60 595454 1 60
61 590865 1 61
62 589379 1 62
63 584428 1 63
64 573100 1 64
65 567456 1 65
66 569028 1 66
67 620735 1 67
68 628884 1 68
69 628232 1 69
70 612117 1 70
71 595404 1 71
72 597141 1 72
73 593408 1 73
74 590072 1 74
75 579799 1 75
76 574205 1 76
77 572775 1 77
78 572942 1 78
79 619567 1 79
80 625809 1 80
81 619916 1 81
82 587625 1 82
83 565742 1 83
84 557274 1 84
85 560576 1 85
86 548854 0 86
87 531673 0 87
88 525919 0 88
89 511038 0 89
90 498662 0 90
91 555362 0 91
92 564591 0 92
93 541657 0 93
94 527070 0 94
95 509846 0 95
96 514258 0 96
97 516922 0 97
98 507561 0 98
99 492622 0 99
100 490243 0 100
101 469357 0 101
102 477580 0 102
103 528379 0 103
104 533590 0 104
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) D t
481854.2 76260.8 503.6
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-63358.6 -22954.2 -104.0 21175.7 75702.4
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 481854.2 6114.0 78.812 < 2e-16 ***
D 76260.8 7172.6 10.632 < 2e-16 ***
t 503.6 109.3 4.608 1.19e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 30950 on 101 degrees of freedom
Multiple R-squared: 0.6649, Adjusted R-squared: 0.6583
F-statistic: 100.2 on 2 and 101 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.006498159 0.01299632 0.99350184
[2,] 0.341404253 0.68280851 0.65859575
[3,] 0.634437785 0.73112443 0.36556221
[4,] 0.539541565 0.92091687 0.46045844
[5,] 0.421866283 0.84373257 0.57813372
[6,] 0.376540337 0.75308067 0.62345966
[7,] 0.328681333 0.65736267 0.67131867
[8,] 0.253836602 0.50767320 0.74616340
[9,] 0.201935448 0.40387090 0.79806455
[10,] 0.183962475 0.36792495 0.81603753
[11,] 0.173142435 0.34628487 0.82685757
[12,] 0.177437547 0.35487509 0.82256245
[13,] 0.194497581 0.38899516 0.80550242
[14,] 0.231893683 0.46378737 0.76810632
[15,] 0.379246900 0.75849380 0.62075310
[16,] 0.469025249 0.93805050 0.53097475
[17,] 0.457733985 0.91546797 0.54226602
[18,] 0.418259433 0.83651887 0.58174057
[19,] 0.385115414 0.77023083 0.61488459
[20,] 0.350167486 0.70033497 0.64983251
[21,] 0.324830986 0.64966197 0.67516901
[22,] 0.317364360 0.63472872 0.68263564
[23,] 0.355553368 0.71110674 0.64444663
[24,] 0.456050226 0.91210045 0.54394977
[25,] 0.595195710 0.80960858 0.40480429
[26,] 0.689784777 0.62043045 0.31021522
[27,] 0.767348259 0.46530348 0.23265174
[28,] 0.799366524 0.40126695 0.20063348
[29,] 0.806406735 0.38718653 0.19359326
[30,] 0.804655896 0.39068821 0.19534410
[31,] 0.802541108 0.39491778 0.19745889
[32,] 0.798802989 0.40239402 0.20119701
[33,] 0.794604640 0.41079072 0.20539536
[34,] 0.803622213 0.39275557 0.19637779
[35,] 0.823917855 0.35216429 0.17608214
[36,] 0.880568040 0.23886392 0.11943196
[37,] 0.925449181 0.14910164 0.07455082
[38,] 0.960521960 0.07895608 0.03947804
[39,] 0.982553690 0.03489262 0.01744631
[40,] 0.989827663 0.02034467 0.01017234
[41,] 0.989548482 0.02090304 0.01045152
[42,] 0.985740798 0.02851840 0.01425920
[43,] 0.980955042 0.03808992 0.01904496
[44,] 0.976385547 0.04722891 0.02361445
[45,] 0.969885928 0.06022814 0.03011407
[46,] 0.959667338 0.08066532 0.04033266
[47,] 0.946111941 0.10777612 0.05388806
[48,] 0.934160184 0.13167963 0.06583982
[49,] 0.917312552 0.16537490 0.08268745
[50,] 0.898626490 0.20274702 0.10137351
[51,] 0.873000477 0.25399905 0.12699952
[52,] 0.842545885 0.31490823 0.15745412
[53,] 0.806919129 0.38616174 0.19308087
[54,] 0.775531778 0.44893644 0.22446822
[55,] 0.738198776 0.52360245 0.26180122
[56,] 0.707140403 0.58571919 0.29285960
[57,] 0.678625986 0.64274803 0.32137401
[58,] 0.667271289 0.66545742 0.33272871
[59,] 0.716597672 0.56680466 0.28340233
[60,] 0.807474610 0.38505078 0.19252539
[61,] 0.887702857 0.22459429 0.11229714
[62,] 0.861275774 0.27744845 0.13872423
[63,] 0.846454362 0.30709128 0.15354564
[64,] 0.835195634 0.32960873 0.16480437
[65,] 0.797528336 0.40494333 0.20247166
[66,] 0.760568212 0.47886358 0.23943179
[67,] 0.716846913 0.56630617 0.28315309
[68,] 0.674852314 0.65029537 0.32514769
[69,] 0.636433075 0.72713385 0.36356693
[70,] 0.630312164 0.73937567 0.36968784
[71,] 0.650779351 0.69844130 0.34922065
[72,] 0.679592605 0.64081479 0.32040740
[73,] 0.709329878 0.58134024 0.29067012
[74,] 0.686946903 0.62610619 0.31305310
[75,] 0.734432818 0.53113436 0.26556718
[76,] 0.813401149 0.37319770 0.18659885
[77,] 0.804565265 0.39086947 0.19543473
[78,] 0.790120324 0.41975935 0.20987968
[79,] 0.781893737 0.43621253 0.21810626
[80,] 0.760149716 0.47970057 0.23985028
[81,] 0.722061593 0.55587681 0.27793841
[82,] 0.677099845 0.64580031 0.32290016
[83,] 0.628989564 0.74202087 0.37101044
[84,] 0.635114771 0.72977046 0.36488523
[85,] 0.760956616 0.47808677 0.23904338
[86,] 0.701507007 0.59698599 0.29849299
[87,] 0.752138764 0.49572247 0.24786124
[88,] 0.736571971 0.52685606 0.26342803
[89,] 0.690728852 0.61854230 0.30927115
[90,] 0.596469669 0.80706066 0.40353033
[91,] 0.514400722 0.97119856 0.48559928
[92,] 0.502840744 0.99431851 0.49715926
[93,] 0.521787316 0.95642537 0.47821268
> postscript(file="/var/www/html/freestat/rcomp/tmp/1lc451227789086.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/2vogp1227789086.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/36htg1227789086.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/43dhs1227789086.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/5vrsq1227789086.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
10507.20604 -1900.37184 -21429.94972 -27260.52760 -42395.10548 -45727.68336
7 8 9 10 11 12
2800.73876 34681.16088 15105.58300 -1864.99488 -23197.57276 -27727.15064
13 14 15 16 17 18
-21363.72852 -28834.30640 -41419.88428 -47044.46216 -54328.04004 -59590.61792
19 20 21 22 23 24
-7407.19580 17747.22632 20497.64844 9898.07056 -22452.50732 -22873.08520
25 26 27 28 29 30
-18394.66309 -20342.24097 -25011.81885 -34703.39673 -41733.97461 -41335.55249
31 32 33 34 35 36
19381.86963 27223.29175 24502.71387 19609.13599 9759.55811 12254.98023
37 38 39 40 41 42
18677.40235 16018.82447 8439.24659 7129.66871 -1643.90917 3966.51295
43 44 45 46 47 48
65814.93507 75702.35719 73476.77931 60445.20143 41821.62355 48762.04567
49 50 51 52 53 54
55795.46779 53820.88990 47795.31202 35558.73414 28118.15626 33674.57838
55 56 57 58 59 60
7718.18433 24447.60645 25794.02857 24001.45069 6340.87281 7124.29493
61 62 63 64 65 66
2031.71705 42.13917 -5412.43871 -17244.01659 -23391.59447 -22323.17235
67 68 69 70 71 72
28880.24977 36525.67189 35370.09401 18751.51613 1534.93825 2768.36037
73 74 75 76 77 78
-1468.21751 -5307.79539 -16084.37327 -22181.95115 -24115.52903 -24452.10691
79 80 81 82 83 84
21669.31521 27407.73733 21011.15944 -11783.41844 -34169.99632 -43141.57420
85 86 87 88 89 90
-40343.15208 23692.08621 6007.50833 -250.06955 -15634.64743 -28514.22531
91 92 93 94 95 96
27682.19681 36407.61893 12970.04105 -2120.53683 -19848.11471 -15939.69259
97 98 99 100 101 102
-13779.27047 -23643.84835 -39086.42623 -41969.00411 -63358.58199 -55639.15988
103 104
-5343.73776 -636.31564
> postscript(file="/var/www/html/freestat/rcomp/tmp/6ab3r1227789086.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 10507.20604 NA
1 -1900.37184 10507.20604
2 -21429.94972 -1900.37184
3 -27260.52760 -21429.94972
4 -42395.10548 -27260.52760
5 -45727.68336 -42395.10548
6 2800.73876 -45727.68336
7 34681.16088 2800.73876
8 15105.58300 34681.16088
9 -1864.99488 15105.58300
10 -23197.57276 -1864.99488
11 -27727.15064 -23197.57276
12 -21363.72852 -27727.15064
13 -28834.30640 -21363.72852
14 -41419.88428 -28834.30640
15 -47044.46216 -41419.88428
16 -54328.04004 -47044.46216
17 -59590.61792 -54328.04004
18 -7407.19580 -59590.61792
19 17747.22632 -7407.19580
20 20497.64844 17747.22632
21 9898.07056 20497.64844
22 -22452.50732 9898.07056
23 -22873.08520 -22452.50732
24 -18394.66309 -22873.08520
25 -20342.24097 -18394.66309
26 -25011.81885 -20342.24097
27 -34703.39673 -25011.81885
28 -41733.97461 -34703.39673
29 -41335.55249 -41733.97461
30 19381.86963 -41335.55249
31 27223.29175 19381.86963
32 24502.71387 27223.29175
33 19609.13599 24502.71387
34 9759.55811 19609.13599
35 12254.98023 9759.55811
36 18677.40235 12254.98023
37 16018.82447 18677.40235
38 8439.24659 16018.82447
39 7129.66871 8439.24659
40 -1643.90917 7129.66871
41 3966.51295 -1643.90917
42 65814.93507 3966.51295
43 75702.35719 65814.93507
44 73476.77931 75702.35719
45 60445.20143 73476.77931
46 41821.62355 60445.20143
47 48762.04567 41821.62355
48 55795.46779 48762.04567
49 53820.88990 55795.46779
50 47795.31202 53820.88990
51 35558.73414 47795.31202
52 28118.15626 35558.73414
53 33674.57838 28118.15626
54 7718.18433 33674.57838
55 24447.60645 7718.18433
56 25794.02857 24447.60645
57 24001.45069 25794.02857
58 6340.87281 24001.45069
59 7124.29493 6340.87281
60 2031.71705 7124.29493
61 42.13917 2031.71705
62 -5412.43871 42.13917
63 -17244.01659 -5412.43871
64 -23391.59447 -17244.01659
65 -22323.17235 -23391.59447
66 28880.24977 -22323.17235
67 36525.67189 28880.24977
68 35370.09401 36525.67189
69 18751.51613 35370.09401
70 1534.93825 18751.51613
71 2768.36037 1534.93825
72 -1468.21751 2768.36037
73 -5307.79539 -1468.21751
74 -16084.37327 -5307.79539
75 -22181.95115 -16084.37327
76 -24115.52903 -22181.95115
77 -24452.10691 -24115.52903
78 21669.31521 -24452.10691
79 27407.73733 21669.31521
80 21011.15944 27407.73733
81 -11783.41844 21011.15944
82 -34169.99632 -11783.41844
83 -43141.57420 -34169.99632
84 -40343.15208 -43141.57420
85 23692.08621 -40343.15208
86 6007.50833 23692.08621
87 -250.06955 6007.50833
88 -15634.64743 -250.06955
89 -28514.22531 -15634.64743
90 27682.19681 -28514.22531
91 36407.61893 27682.19681
92 12970.04105 36407.61893
93 -2120.53683 12970.04105
94 -19848.11471 -2120.53683
95 -15939.69259 -19848.11471
96 -13779.27047 -15939.69259
97 -23643.84835 -13779.27047
98 -39086.42623 -23643.84835
99 -41969.00411 -39086.42623
100 -63358.58199 -41969.00411
101 -55639.15988 -63358.58199
102 -5343.73776 -55639.15988
103 -636.31564 -5343.73776
104 NA -636.31564
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1900.37184 10507.20604
[2,] -21429.94972 -1900.37184
[3,] -27260.52760 -21429.94972
[4,] -42395.10548 -27260.52760
[5,] -45727.68336 -42395.10548
[6,] 2800.73876 -45727.68336
[7,] 34681.16088 2800.73876
[8,] 15105.58300 34681.16088
[9,] -1864.99488 15105.58300
[10,] -23197.57276 -1864.99488
[11,] -27727.15064 -23197.57276
[12,] -21363.72852 -27727.15064
[13,] -28834.30640 -21363.72852
[14,] -41419.88428 -28834.30640
[15,] -47044.46216 -41419.88428
[16,] -54328.04004 -47044.46216
[17,] -59590.61792 -54328.04004
[18,] -7407.19580 -59590.61792
[19,] 17747.22632 -7407.19580
[20,] 20497.64844 17747.22632
[21,] 9898.07056 20497.64844
[22,] -22452.50732 9898.07056
[23,] -22873.08520 -22452.50732
[24,] -18394.66309 -22873.08520
[25,] -20342.24097 -18394.66309
[26,] -25011.81885 -20342.24097
[27,] -34703.39673 -25011.81885
[28,] -41733.97461 -34703.39673
[29,] -41335.55249 -41733.97461
[30,] 19381.86963 -41335.55249
[31,] 27223.29175 19381.86963
[32,] 24502.71387 27223.29175
[33,] 19609.13599 24502.71387
[34,] 9759.55811 19609.13599
[35,] 12254.98023 9759.55811
[36,] 18677.40235 12254.98023
[37,] 16018.82447 18677.40235
[38,] 8439.24659 16018.82447
[39,] 7129.66871 8439.24659
[40,] -1643.90917 7129.66871
[41,] 3966.51295 -1643.90917
[42,] 65814.93507 3966.51295
[43,] 75702.35719 65814.93507
[44,] 73476.77931 75702.35719
[45,] 60445.20143 73476.77931
[46,] 41821.62355 60445.20143
[47,] 48762.04567 41821.62355
[48,] 55795.46779 48762.04567
[49,] 53820.88990 55795.46779
[50,] 47795.31202 53820.88990
[51,] 35558.73414 47795.31202
[52,] 28118.15626 35558.73414
[53,] 33674.57838 28118.15626
[54,] 7718.18433 33674.57838
[55,] 24447.60645 7718.18433
[56,] 25794.02857 24447.60645
[57,] 24001.45069 25794.02857
[58,] 6340.87281 24001.45069
[59,] 7124.29493 6340.87281
[60,] 2031.71705 7124.29493
[61,] 42.13917 2031.71705
[62,] -5412.43871 42.13917
[63,] -17244.01659 -5412.43871
[64,] -23391.59447 -17244.01659
[65,] -22323.17235 -23391.59447
[66,] 28880.24977 -22323.17235
[67,] 36525.67189 28880.24977
[68,] 35370.09401 36525.67189
[69,] 18751.51613 35370.09401
[70,] 1534.93825 18751.51613
[71,] 2768.36037 1534.93825
[72,] -1468.21751 2768.36037
[73,] -5307.79539 -1468.21751
[74,] -16084.37327 -5307.79539
[75,] -22181.95115 -16084.37327
[76,] -24115.52903 -22181.95115
[77,] -24452.10691 -24115.52903
[78,] 21669.31521 -24452.10691
[79,] 27407.73733 21669.31521
[80,] 21011.15944 27407.73733
[81,] -11783.41844 21011.15944
[82,] -34169.99632 -11783.41844
[83,] -43141.57420 -34169.99632
[84,] -40343.15208 -43141.57420
[85,] 23692.08621 -40343.15208
[86,] 6007.50833 23692.08621
[87,] -250.06955 6007.50833
[88,] -15634.64743 -250.06955
[89,] -28514.22531 -15634.64743
[90,] 27682.19681 -28514.22531
[91,] 36407.61893 27682.19681
[92,] 12970.04105 36407.61893
[93,] -2120.53683 12970.04105
[94,] -19848.11471 -2120.53683
[95,] -15939.69259 -19848.11471
[96,] -13779.27047 -15939.69259
[97,] -23643.84835 -13779.27047
[98,] -39086.42623 -23643.84835
[99,] -41969.00411 -39086.42623
[100,] -63358.58199 -41969.00411
[101,] -55639.15988 -63358.58199
[102,] -5343.73776 -55639.15988
[103,] -636.31564 -5343.73776
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1900.37184 10507.20604
2 -21429.94972 -1900.37184
3 -27260.52760 -21429.94972
4 -42395.10548 -27260.52760
5 -45727.68336 -42395.10548
6 2800.73876 -45727.68336
7 34681.16088 2800.73876
8 15105.58300 34681.16088
9 -1864.99488 15105.58300
10 -23197.57276 -1864.99488
11 -27727.15064 -23197.57276
12 -21363.72852 -27727.15064
13 -28834.30640 -21363.72852
14 -41419.88428 -28834.30640
15 -47044.46216 -41419.88428
16 -54328.04004 -47044.46216
17 -59590.61792 -54328.04004
18 -7407.19580 -59590.61792
19 17747.22632 -7407.19580
20 20497.64844 17747.22632
21 9898.07056 20497.64844
22 -22452.50732 9898.07056
23 -22873.08520 -22452.50732
24 -18394.66309 -22873.08520
25 -20342.24097 -18394.66309
26 -25011.81885 -20342.24097
27 -34703.39673 -25011.81885
28 -41733.97461 -34703.39673
29 -41335.55249 -41733.97461
30 19381.86963 -41335.55249
31 27223.29175 19381.86963
32 24502.71387 27223.29175
33 19609.13599 24502.71387
34 9759.55811 19609.13599
35 12254.98023 9759.55811
36 18677.40235 12254.98023
37 16018.82447 18677.40235
38 8439.24659 16018.82447
39 7129.66871 8439.24659
40 -1643.90917 7129.66871
41 3966.51295 -1643.90917
42 65814.93507 3966.51295
43 75702.35719 65814.93507
44 73476.77931 75702.35719
45 60445.20143 73476.77931
46 41821.62355 60445.20143
47 48762.04567 41821.62355
48 55795.46779 48762.04567
49 53820.88990 55795.46779
50 47795.31202 53820.88990
51 35558.73414 47795.31202
52 28118.15626 35558.73414
53 33674.57838 28118.15626
54 7718.18433 33674.57838
55 24447.60645 7718.18433
56 25794.02857 24447.60645
57 24001.45069 25794.02857
58 6340.87281 24001.45069
59 7124.29493 6340.87281
60 2031.71705 7124.29493
61 42.13917 2031.71705
62 -5412.43871 42.13917
63 -17244.01659 -5412.43871
64 -23391.59447 -17244.01659
65 -22323.17235 -23391.59447
66 28880.24977 -22323.17235
67 36525.67189 28880.24977
68 35370.09401 36525.67189
69 18751.51613 35370.09401
70 1534.93825 18751.51613
71 2768.36037 1534.93825
72 -1468.21751 2768.36037
73 -5307.79539 -1468.21751
74 -16084.37327 -5307.79539
75 -22181.95115 -16084.37327
76 -24115.52903 -22181.95115
77 -24452.10691 -24115.52903
78 21669.31521 -24452.10691
79 27407.73733 21669.31521
80 21011.15944 27407.73733
81 -11783.41844 21011.15944
82 -34169.99632 -11783.41844
83 -43141.57420 -34169.99632
84 -40343.15208 -43141.57420
85 23692.08621 -40343.15208
86 6007.50833 23692.08621
87 -250.06955 6007.50833
88 -15634.64743 -250.06955
89 -28514.22531 -15634.64743
90 27682.19681 -28514.22531
91 36407.61893 27682.19681
92 12970.04105 36407.61893
93 -2120.53683 12970.04105
94 -19848.11471 -2120.53683
95 -15939.69259 -19848.11471
96 -13779.27047 -15939.69259
97 -23643.84835 -13779.27047
98 -39086.42623 -23643.84835
99 -41969.00411 -39086.42623
100 -63358.58199 -41969.00411
101 -55639.15988 -63358.58199
102 -5343.73776 -55639.15988
103 -636.31564 -5343.73776
> 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/7d2ro1227789086.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/8fyw81227789086.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/9pgqu1227789086.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/10vb6v1227789086.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/11s6161227789086.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/12u4qj1227789087.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/13k1dp1227789087.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/143ytu1227789087.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/15llk01227789087.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/16c5vj1227789087.tab")
+ }
>
> system("convert tmp/1lc451227789086.ps tmp/1lc451227789086.png")
> system("convert tmp/2vogp1227789086.ps tmp/2vogp1227789086.png")
> system("convert tmp/36htg1227789086.ps tmp/36htg1227789086.png")
> system("convert tmp/43dhs1227789086.ps tmp/43dhs1227789086.png")
> system("convert tmp/5vrsq1227789086.ps tmp/5vrsq1227789086.png")
> system("convert tmp/6ab3r1227789086.ps tmp/6ab3r1227789086.png")
> system("convert tmp/7d2ro1227789086.ps tmp/7d2ro1227789086.png")
> system("convert tmp/8fyw81227789086.ps tmp/8fyw81227789086.png")
> system("convert tmp/9pgqu1227789086.ps tmp/9pgqu1227789086.png")
> system("convert tmp/10vb6v1227789086.ps tmp/10vb6v1227789086.png")
>
>
> proc.time()
user system elapsed
4.326 2.638 4.742