R version 2.7.0 (2008-04-22)
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 = '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
W D M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 492865 0 1 0 0 0 0 0 0 0 0 0 0 1
2 480961 0 0 1 0 0 0 0 0 0 0 0 0 2
3 461935 0 0 0 1 0 0 0 0 0 0 0 0 3
4 456608 0 0 0 0 1 0 0 0 0 0 0 0 4
5 441977 0 0 0 0 0 1 0 0 0 0 0 0 5
6 439148 0 0 0 0 0 0 1 0 0 0 0 0 6
7 488180 0 0 0 0 0 0 0 1 0 0 0 0 7
8 520564 0 0 0 0 0 0 0 0 1 0 0 0 8
9 501492 0 0 0 0 0 0 0 0 0 1 0 0 9
10 485025 0 0 0 0 0 0 0 0 0 0 1 0 10
11 464196 0 0 0 0 0 0 0 0 0 0 0 1 11
12 460170 0 0 0 0 0 0 0 0 0 0 0 0 12
13 467037 0 1 0 0 0 0 0 0 0 0 0 0 13
14 460070 0 0 1 0 0 0 0 0 0 0 0 0 14
15 447988 0 0 0 1 0 0 0 0 0 0 0 0 15
16 442867 0 0 0 0 1 0 0 0 0 0 0 0 16
17 436087 0 0 0 0 0 1 0 0 0 0 0 0 17
18 431328 0 0 0 0 0 0 1 0 0 0 0 0 18
19 484015 0 0 0 0 0 0 0 1 0 0 0 0 19
20 509673 0 0 0 0 0 0 0 0 1 0 0 0 20
21 512927 0 0 0 0 0 0 0 0 0 1 0 0 21
22 502831 0 0 0 0 0 0 0 0 0 0 1 0 22
23 470984 0 0 0 0 0 0 0 0 0 0 0 1 23
24 471067 0 0 0 0 0 0 0 0 0 0 0 0 24
25 476049 0 1 0 0 0 0 0 0 0 0 0 0 25
26 474605 0 0 1 0 0 0 0 0 0 0 0 0 26
27 470439 0 0 0 1 0 0 0 0 0 0 0 0 27
28 461251 0 0 0 0 1 0 0 0 0 0 0 0 28
29 454724 0 0 0 0 0 1 0 0 0 0 0 0 29
30 455626 0 0 0 0 0 0 1 0 0 0 0 0 30
31 516847 0 0 0 0 0 0 0 1 0 0 0 0 31
32 525192 0 0 0 0 0 0 0 0 1 0 0 0 32
33 522975 0 0 0 0 0 0 0 0 0 1 0 0 33
34 518585 0 0 0 0 0 0 0 0 0 0 1 0 34
35 509239 0 0 0 0 0 0 0 0 0 0 0 1 35
36 512238 0 0 0 0 0 0 0 0 0 0 0 0 36
37 519164 0 1 0 0 0 0 0 0 0 0 0 0 37
38 517009 0 0 1 0 0 0 0 0 0 0 0 0 38
39 509933 0 0 0 1 0 0 0 0 0 0 0 0 39
40 509127 0 0 0 0 1 0 0 0 0 0 0 0 40
41 500857 0 0 0 0 0 1 0 0 0 0 0 0 41
42 506971 0 0 0 0 0 0 1 0 0 0 0 0 42
43 569323 0 0 0 0 0 0 0 1 0 0 0 0 43
44 579714 0 0 0 0 0 0 0 0 1 0 0 0 44
45 577992 0 0 0 0 0 0 0 0 0 1 0 0 45
46 565464 0 0 0 0 0 0 0 0 0 0 1 0 46
47 547344 0 0 0 0 0 0 0 0 0 0 0 1 47
48 554788 0 0 0 0 0 0 0 0 0 0 0 0 48
49 562325 0 1 0 0 0 0 0 0 0 0 0 0 49
50 560854 0 0 1 0 0 0 0 0 0 0 0 0 50
51 555332 0 0 0 1 0 0 0 0 0 0 0 0 51
52 543599 0 0 0 0 1 0 0 0 0 0 0 0 52
53 536662 0 0 0 0 0 1 0 0 0 0 0 0 53
54 542722 1 0 0 0 0 0 1 0 0 0 0 0 54
55 593530 1 0 0 0 0 0 0 1 0 0 0 0 55
56 610763 1 0 0 0 0 0 0 0 1 0 0 0 56
57 612613 1 0 0 0 0 0 0 0 0 1 0 0 57
58 611324 1 0 0 0 0 0 0 0 0 0 1 0 58
59 594167 1 0 0 0 0 0 0 0 0 0 0 1 59
60 595454 1 0 0 0 0 0 0 0 0 0 0 0 60
61 590865 1 1 0 0 0 0 0 0 0 0 0 0 61
62 589379 1 0 1 0 0 0 0 0 0 0 0 0 62
63 584428 1 0 0 1 0 0 0 0 0 0 0 0 63
64 573100 1 0 0 0 1 0 0 0 0 0 0 0 64
65 567456 1 0 0 0 0 1 0 0 0 0 0 0 65
66 569028 1 0 0 0 0 0 1 0 0 0 0 0 66
67 620735 1 0 0 0 0 0 0 1 0 0 0 0 67
68 628884 1 0 0 0 0 0 0 0 1 0 0 0 68
69 628232 1 0 0 0 0 0 0 0 0 1 0 0 69
70 612117 1 0 0 0 0 0 0 0 0 0 1 0 70
71 595404 1 0 0 0 0 0 0 0 0 0 0 1 71
72 597141 1 0 0 0 0 0 0 0 0 0 0 0 72
73 593408 1 1 0 0 0 0 0 0 0 0 0 0 73
74 590072 1 0 1 0 0 0 0 0 0 0 0 0 74
75 579799 1 0 0 1 0 0 0 0 0 0 0 0 75
76 574205 1 0 0 0 1 0 0 0 0 0 0 0 76
77 572775 1 0 0 0 0 1 0 0 0 0 0 0 77
78 572942 1 0 0 0 0 0 1 0 0 0 0 0 78
79 619567 1 0 0 0 0 0 0 1 0 0 0 0 79
80 625809 1 0 0 0 0 0 0 0 1 0 0 0 80
81 619916 1 0 0 0 0 0 0 0 0 1 0 0 81
82 587625 1 0 0 0 0 0 0 0 0 0 1 0 82
83 565742 1 0 0 0 0 0 0 0 0 0 0 1 83
84 557274 1 0 0 0 0 0 0 0 0 0 0 0 84
85 560576 1 1 0 0 0 0 0 0 0 0 0 0 85
86 548854 0 0 1 0 0 0 0 0 0 0 0 0 86
87 531673 0 0 0 1 0 0 0 0 0 0 0 0 87
88 525919 0 0 0 0 1 0 0 0 0 0 0 0 88
89 511038 0 0 0 0 0 1 0 0 0 0 0 0 89
90 498662 0 0 0 0 0 0 1 0 0 0 0 0 90
91 555362 0 0 0 0 0 0 0 1 0 0 0 0 91
92 564591 0 0 0 0 0 0 0 0 1 0 0 0 92
93 541657 0 0 0 0 0 0 0 0 0 1 0 0 93
94 527070 0 0 0 0 0 0 0 0 0 0 1 0 94
95 509846 0 0 0 0 0 0 0 0 0 0 0 1 95
96 514258 0 0 0 0 0 0 0 0 0 0 0 0 96
97 516922 0 1 0 0 0 0 0 0 0 0 0 0 97
98 507561 0 0 1 0 0 0 0 0 0 0 0 0 98
99 492622 0 0 0 1 0 0 0 0 0 0 0 0 99
100 490243 0 0 0 0 1 0 0 0 0 0 0 0 100
101 469357 0 0 0 0 0 1 0 0 0 0 0 0 101
102 477580 0 0 0 0 0 0 1 0 0 0 0 0 102
103 528379 0 0 0 0 0 0 0 1 0 0 0 0 103
104 533590 0 0 0 0 0 0 0 0 1 0 0 0 104
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) D M1 M2 M3 M4
477621.4 72190.5 3835.0 5797.3 -5302.8 -12182.1
M5 M6 M7 M8 M9 M10
-22256.6 -30456.7 22570.7 35699.3 33488.2 19497.3
M11 t
-163.0 520.5
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-38576 -17998 -3240 13112 56469
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 477621.44 9883.83 48.324 < 2e-16 ***
D 72190.51 5794.28 12.459 < 2e-16 ***
M1 3835.03 12148.88 0.316 0.75298
M2 5797.27 12169.38 0.476 0.63496
M3 -5302.76 12169.57 -0.436 0.66407
M4 -12182.13 12170.40 -1.001 0.31953
M5 -22256.61 12171.88 -1.829 0.07078 .
M6 -30456.70 12144.17 -2.508 0.01394 *
M7 22570.71 12145.15 1.858 0.06638 .
M8 35699.34 12146.79 2.939 0.00418 **
M9 33488.19 12496.59 2.680 0.00876 **
M10 19497.33 12495.02 1.560 0.12217
M11 -163.02 12494.09 -0.013 0.98962
t 520.48 88.41 5.887 6.65e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 24990 on 90 degrees of freedom
Multiple R-squared: 0.8054, Adjusted R-squared: 0.7773
F-statistic: 28.65 on 13 and 90 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.0207481872 4.149637e-02 9.792518e-01
[2,] 0.0068123226 1.362465e-02 9.931877e-01
[3,] 0.0030767543 6.153509e-03 9.969232e-01
[4,] 0.0007494608 1.498922e-03 9.992505e-01
[5,] 0.0027426571 5.485314e-03 9.972573e-01
[6,] 0.0065551344 1.311027e-02 9.934449e-01
[7,] 0.0043338477 8.667695e-03 9.956662e-01
[8,] 0.0036233969 7.246794e-03 9.963766e-01
[9,] 0.0018495549 3.699110e-03 9.981504e-01
[10,] 0.0014154672 2.830934e-03 9.985845e-01
[11,] 0.0022626730 4.525346e-03 9.977373e-01
[12,] 0.0025071880 5.014376e-03 9.974928e-01
[13,] 0.0035792885 7.158577e-03 9.964207e-01
[14,] 0.0057887942 1.157759e-02 9.942112e-01
[15,] 0.0157758856 3.155177e-02 9.842241e-01
[16,] 0.0155119894 3.102398e-02 9.844880e-01
[17,] 0.0192698582 3.853972e-02 9.807301e-01
[18,] 0.0286550006 5.731000e-02 9.713450e-01
[19,] 0.0868547094 1.737094e-01 9.131453e-01
[20,] 0.2128193729 4.256387e-01 7.871806e-01
[21,] 0.2741001381 5.482003e-01 7.258999e-01
[22,] 0.4113962165 8.227924e-01 5.886038e-01
[23,] 0.5823690866 8.352618e-01 4.176309e-01
[24,] 0.7440510377 5.118979e-01 2.559490e-01
[25,] 0.8639975011 2.720050e-01 1.360025e-01
[26,] 0.9245065636 1.509869e-01 7.549344e-02
[27,] 0.9600986159 7.980277e-02 3.990138e-02
[28,] 0.9635512251 7.289755e-02 3.644877e-02
[29,] 0.9677819903 6.443602e-02 3.221801e-02
[30,] 0.9668431087 6.631378e-02 3.315689e-02
[31,] 0.9666754008 6.664920e-02 3.332460e-02
[32,] 0.9699307782 6.013844e-02 3.006922e-02
[33,] 0.9654895945 6.902081e-02 3.451041e-02
[34,] 0.9603329085 7.933418e-02 3.966709e-02
[35,] 0.9596758332 8.064833e-02 4.032417e-02
[36,] 0.9507012744 9.859745e-02 4.929873e-02
[37,] 0.9459102705 1.081795e-01 5.408973e-02
[38,] 0.9678826098 6.423478e-02 3.211739e-02
[39,] 0.9890268266 2.194635e-02 1.097317e-02
[40,] 0.9957766411 8.446718e-03 4.223359e-03
[41,] 0.9979006994 4.198601e-03 2.099301e-03
[42,] 0.9971423281 5.715344e-03 2.857672e-03
[43,] 0.9961058401 7.788320e-03 3.894160e-03
[44,] 0.9944568967 1.108621e-02 5.543103e-03
[45,] 0.9940895951 1.182081e-02 5.910405e-03
[46,] 0.9963087492 7.382502e-03 3.691251e-03
[47,] 0.9961128318 7.774336e-03 3.887168e-03
[48,] 0.9981997556 3.600489e-03 1.800244e-03
[49,] 0.9988241562 2.351688e-03 1.175844e-03
[50,] 0.9992548559 1.490288e-03 7.451441e-04
[51,] 0.9997103927 5.792145e-04 2.896073e-04
[52,] 0.9999560254 8.794925e-05 4.397463e-05
[53,] 0.9999557759 8.844816e-05 4.422408e-05
[54,] 0.9999206886 1.586228e-04 7.931139e-05
[55,] 0.9998341233 3.317535e-04 1.658767e-04
[56,] 0.9996350502 7.298996e-04 3.649498e-04
[57,] 0.9994983666 1.003267e-03 5.016334e-04
[58,] 0.9995612547 8.774906e-04 4.387453e-04
[59,] 0.9993717143 1.256571e-03 6.282857e-04
[60,] 0.9992613211 1.477358e-03 7.386789e-04
[61,] 0.9985917904 2.816419e-03 1.408210e-03
[62,] 0.9979212492 4.157502e-03 2.078751e-03
[63,] 0.9957453018 8.509396e-03 4.254698e-03
[64,] 0.9921530890 1.569382e-02 7.846911e-03
[65,] 0.9989479630 2.104074e-03 1.052037e-03
[66,] 0.9991246048 1.750790e-03 8.753952e-04
[67,] 0.9990892943 1.821411e-03 9.107057e-04
[68,] 0.9976499722 4.700056e-03 2.350028e-03
[69,] 0.9935881560 1.282369e-02 6.411844e-03
[70,] 0.9856781723 2.864366e-02 1.432183e-02
[71,] 0.9626497043 7.470059e-02 3.735030e-02
> postscript(file="/var/www/html/rcomp/tmp/18os01227791115.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/229e21227791115.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/3iima1227791115.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/4kcny1227791115.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/5fme81227791115.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
10888.0501 -3498.6736 -11945.1180 -10913.2291 -15990.2291 -11139.6165
7 8 9 10 11 12
-15655.5054 3079.3835 -14301.9411 -17298.5661 -18987.6911 -23697.1911
13 14 15 16 17 18
-21185.6974 -30635.4211 -32137.8656 -30899.9767 -28125.9767 -25205.3641
19 20 21 22 23 24
-26066.2529 -14057.3641 -9112.6887 -5738.3137 -18445.4387 -19045.9387
25 26 27 28 29 30
-18419.4449 -22346.1686 -15932.6131 -18761.7242 -15734.7242 -7153.1116
31 32 33 34 35 36
519.9995 -4784.1116 -5310.4362 3769.9388 13563.8138 15879.3138
37 38 39 40 41 42
18449.8076 13812.0838 17315.6394 22868.5283 24152.5283 37946.1409
43 44 45 46 47 48
46750.2520 43492.1409 43460.8163 44403.1913 45423.0663 52183.5663
49 50 51 52 53 54
55365.0600 51411.3363 56468.8919 51094.7808 53711.7808 -4739.1201
55 56 57 58 59 60
-7479.0090 -3895.1201 -354.4447 11826.9303 13809.8053 14413.3053
61 62 63 64 65 66
5468.7991 1500.0754 7128.6309 2159.5198 6069.5198 15321.1324
67 68 69 70 71 72
13480.2435 7980.1324 9018.8078 6374.1828 8801.0578 9854.5578
73 74 75 76 77 78
1766.0515 -4052.6722 -3746.1166 -2981.2277 5142.7723 12989.3849
79 80 81 82 83 84
6066.4960 -1340.6151 -5542.9397 -24363.5647 -27106.6897 -36258.1897
85 86 87 88 89 90
-37311.6960 20674.0938 14072.6493 14677.5382 9350.5382 4654.1508
91 92 93 94 95 96
7806.2619 3386.1508 -17857.1738 -18973.7988 -17057.9238 -13329.4238
97 98 99 100 101 102
-15020.9301 -26864.6538 -31224.0982 -27244.2093 -38576.2093 -22673.5967
103 104
-25422.4856 -33860.5967
> postscript(file="/var/www/html/rcomp/tmp/6g4b31227791115.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 10888.0501 NA
1 -3498.6736 10888.0501
2 -11945.1180 -3498.6736
3 -10913.2291 -11945.1180
4 -15990.2291 -10913.2291
5 -11139.6165 -15990.2291
6 -15655.5054 -11139.6165
7 3079.3835 -15655.5054
8 -14301.9411 3079.3835
9 -17298.5661 -14301.9411
10 -18987.6911 -17298.5661
11 -23697.1911 -18987.6911
12 -21185.6974 -23697.1911
13 -30635.4211 -21185.6974
14 -32137.8656 -30635.4211
15 -30899.9767 -32137.8656
16 -28125.9767 -30899.9767
17 -25205.3641 -28125.9767
18 -26066.2529 -25205.3641
19 -14057.3641 -26066.2529
20 -9112.6887 -14057.3641
21 -5738.3137 -9112.6887
22 -18445.4387 -5738.3137
23 -19045.9387 -18445.4387
24 -18419.4449 -19045.9387
25 -22346.1686 -18419.4449
26 -15932.6131 -22346.1686
27 -18761.7242 -15932.6131
28 -15734.7242 -18761.7242
29 -7153.1116 -15734.7242
30 519.9995 -7153.1116
31 -4784.1116 519.9995
32 -5310.4362 -4784.1116
33 3769.9388 -5310.4362
34 13563.8138 3769.9388
35 15879.3138 13563.8138
36 18449.8076 15879.3138
37 13812.0838 18449.8076
38 17315.6394 13812.0838
39 22868.5283 17315.6394
40 24152.5283 22868.5283
41 37946.1409 24152.5283
42 46750.2520 37946.1409
43 43492.1409 46750.2520
44 43460.8163 43492.1409
45 44403.1913 43460.8163
46 45423.0663 44403.1913
47 52183.5663 45423.0663
48 55365.0600 52183.5663
49 51411.3363 55365.0600
50 56468.8919 51411.3363
51 51094.7808 56468.8919
52 53711.7808 51094.7808
53 -4739.1201 53711.7808
54 -7479.0090 -4739.1201
55 -3895.1201 -7479.0090
56 -354.4447 -3895.1201
57 11826.9303 -354.4447
58 13809.8053 11826.9303
59 14413.3053 13809.8053
60 5468.7991 14413.3053
61 1500.0754 5468.7991
62 7128.6309 1500.0754
63 2159.5198 7128.6309
64 6069.5198 2159.5198
65 15321.1324 6069.5198
66 13480.2435 15321.1324
67 7980.1324 13480.2435
68 9018.8078 7980.1324
69 6374.1828 9018.8078
70 8801.0578 6374.1828
71 9854.5578 8801.0578
72 1766.0515 9854.5578
73 -4052.6722 1766.0515
74 -3746.1166 -4052.6722
75 -2981.2277 -3746.1166
76 5142.7723 -2981.2277
77 12989.3849 5142.7723
78 6066.4960 12989.3849
79 -1340.6151 6066.4960
80 -5542.9397 -1340.6151
81 -24363.5647 -5542.9397
82 -27106.6897 -24363.5647
83 -36258.1897 -27106.6897
84 -37311.6960 -36258.1897
85 20674.0938 -37311.6960
86 14072.6493 20674.0938
87 14677.5382 14072.6493
88 9350.5382 14677.5382
89 4654.1508 9350.5382
90 7806.2619 4654.1508
91 3386.1508 7806.2619
92 -17857.1738 3386.1508
93 -18973.7988 -17857.1738
94 -17057.9238 -18973.7988
95 -13329.4238 -17057.9238
96 -15020.9301 -13329.4238
97 -26864.6538 -15020.9301
98 -31224.0982 -26864.6538
99 -27244.2093 -31224.0982
100 -38576.2093 -27244.2093
101 -22673.5967 -38576.2093
102 -25422.4856 -22673.5967
103 -33860.5967 -25422.4856
104 NA -33860.5967
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -3498.6736 10888.0501
[2,] -11945.1180 -3498.6736
[3,] -10913.2291 -11945.1180
[4,] -15990.2291 -10913.2291
[5,] -11139.6165 -15990.2291
[6,] -15655.5054 -11139.6165
[7,] 3079.3835 -15655.5054
[8,] -14301.9411 3079.3835
[9,] -17298.5661 -14301.9411
[10,] -18987.6911 -17298.5661
[11,] -23697.1911 -18987.6911
[12,] -21185.6974 -23697.1911
[13,] -30635.4211 -21185.6974
[14,] -32137.8656 -30635.4211
[15,] -30899.9767 -32137.8656
[16,] -28125.9767 -30899.9767
[17,] -25205.3641 -28125.9767
[18,] -26066.2529 -25205.3641
[19,] -14057.3641 -26066.2529
[20,] -9112.6887 -14057.3641
[21,] -5738.3137 -9112.6887
[22,] -18445.4387 -5738.3137
[23,] -19045.9387 -18445.4387
[24,] -18419.4449 -19045.9387
[25,] -22346.1686 -18419.4449
[26,] -15932.6131 -22346.1686
[27,] -18761.7242 -15932.6131
[28,] -15734.7242 -18761.7242
[29,] -7153.1116 -15734.7242
[30,] 519.9995 -7153.1116
[31,] -4784.1116 519.9995
[32,] -5310.4362 -4784.1116
[33,] 3769.9388 -5310.4362
[34,] 13563.8138 3769.9388
[35,] 15879.3138 13563.8138
[36,] 18449.8076 15879.3138
[37,] 13812.0838 18449.8076
[38,] 17315.6394 13812.0838
[39,] 22868.5283 17315.6394
[40,] 24152.5283 22868.5283
[41,] 37946.1409 24152.5283
[42,] 46750.2520 37946.1409
[43,] 43492.1409 46750.2520
[44,] 43460.8163 43492.1409
[45,] 44403.1913 43460.8163
[46,] 45423.0663 44403.1913
[47,] 52183.5663 45423.0663
[48,] 55365.0600 52183.5663
[49,] 51411.3363 55365.0600
[50,] 56468.8919 51411.3363
[51,] 51094.7808 56468.8919
[52,] 53711.7808 51094.7808
[53,] -4739.1201 53711.7808
[54,] -7479.0090 -4739.1201
[55,] -3895.1201 -7479.0090
[56,] -354.4447 -3895.1201
[57,] 11826.9303 -354.4447
[58,] 13809.8053 11826.9303
[59,] 14413.3053 13809.8053
[60,] 5468.7991 14413.3053
[61,] 1500.0754 5468.7991
[62,] 7128.6309 1500.0754
[63,] 2159.5198 7128.6309
[64,] 6069.5198 2159.5198
[65,] 15321.1324 6069.5198
[66,] 13480.2435 15321.1324
[67,] 7980.1324 13480.2435
[68,] 9018.8078 7980.1324
[69,] 6374.1828 9018.8078
[70,] 8801.0578 6374.1828
[71,] 9854.5578 8801.0578
[72,] 1766.0515 9854.5578
[73,] -4052.6722 1766.0515
[74,] -3746.1166 -4052.6722
[75,] -2981.2277 -3746.1166
[76,] 5142.7723 -2981.2277
[77,] 12989.3849 5142.7723
[78,] 6066.4960 12989.3849
[79,] -1340.6151 6066.4960
[80,] -5542.9397 -1340.6151
[81,] -24363.5647 -5542.9397
[82,] -27106.6897 -24363.5647
[83,] -36258.1897 -27106.6897
[84,] -37311.6960 -36258.1897
[85,] 20674.0938 -37311.6960
[86,] 14072.6493 20674.0938
[87,] 14677.5382 14072.6493
[88,] 9350.5382 14677.5382
[89,] 4654.1508 9350.5382
[90,] 7806.2619 4654.1508
[91,] 3386.1508 7806.2619
[92,] -17857.1738 3386.1508
[93,] -18973.7988 -17857.1738
[94,] -17057.9238 -18973.7988
[95,] -13329.4238 -17057.9238
[96,] -15020.9301 -13329.4238
[97,] -26864.6538 -15020.9301
[98,] -31224.0982 -26864.6538
[99,] -27244.2093 -31224.0982
[100,] -38576.2093 -27244.2093
[101,] -22673.5967 -38576.2093
[102,] -25422.4856 -22673.5967
[103,] -33860.5967 -25422.4856
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -3498.6736 10888.0501
2 -11945.1180 -3498.6736
3 -10913.2291 -11945.1180
4 -15990.2291 -10913.2291
5 -11139.6165 -15990.2291
6 -15655.5054 -11139.6165
7 3079.3835 -15655.5054
8 -14301.9411 3079.3835
9 -17298.5661 -14301.9411
10 -18987.6911 -17298.5661
11 -23697.1911 -18987.6911
12 -21185.6974 -23697.1911
13 -30635.4211 -21185.6974
14 -32137.8656 -30635.4211
15 -30899.9767 -32137.8656
16 -28125.9767 -30899.9767
17 -25205.3641 -28125.9767
18 -26066.2529 -25205.3641
19 -14057.3641 -26066.2529
20 -9112.6887 -14057.3641
21 -5738.3137 -9112.6887
22 -18445.4387 -5738.3137
23 -19045.9387 -18445.4387
24 -18419.4449 -19045.9387
25 -22346.1686 -18419.4449
26 -15932.6131 -22346.1686
27 -18761.7242 -15932.6131
28 -15734.7242 -18761.7242
29 -7153.1116 -15734.7242
30 519.9995 -7153.1116
31 -4784.1116 519.9995
32 -5310.4362 -4784.1116
33 3769.9388 -5310.4362
34 13563.8138 3769.9388
35 15879.3138 13563.8138
36 18449.8076 15879.3138
37 13812.0838 18449.8076
38 17315.6394 13812.0838
39 22868.5283 17315.6394
40 24152.5283 22868.5283
41 37946.1409 24152.5283
42 46750.2520 37946.1409
43 43492.1409 46750.2520
44 43460.8163 43492.1409
45 44403.1913 43460.8163
46 45423.0663 44403.1913
47 52183.5663 45423.0663
48 55365.0600 52183.5663
49 51411.3363 55365.0600
50 56468.8919 51411.3363
51 51094.7808 56468.8919
52 53711.7808 51094.7808
53 -4739.1201 53711.7808
54 -7479.0090 -4739.1201
55 -3895.1201 -7479.0090
56 -354.4447 -3895.1201
57 11826.9303 -354.4447
58 13809.8053 11826.9303
59 14413.3053 13809.8053
60 5468.7991 14413.3053
61 1500.0754 5468.7991
62 7128.6309 1500.0754
63 2159.5198 7128.6309
64 6069.5198 2159.5198
65 15321.1324 6069.5198
66 13480.2435 15321.1324
67 7980.1324 13480.2435
68 9018.8078 7980.1324
69 6374.1828 9018.8078
70 8801.0578 6374.1828
71 9854.5578 8801.0578
72 1766.0515 9854.5578
73 -4052.6722 1766.0515
74 -3746.1166 -4052.6722
75 -2981.2277 -3746.1166
76 5142.7723 -2981.2277
77 12989.3849 5142.7723
78 6066.4960 12989.3849
79 -1340.6151 6066.4960
80 -5542.9397 -1340.6151
81 -24363.5647 -5542.9397
82 -27106.6897 -24363.5647
83 -36258.1897 -27106.6897
84 -37311.6960 -36258.1897
85 20674.0938 -37311.6960
86 14072.6493 20674.0938
87 14677.5382 14072.6493
88 9350.5382 14677.5382
89 4654.1508 9350.5382
90 7806.2619 4654.1508
91 3386.1508 7806.2619
92 -17857.1738 3386.1508
93 -18973.7988 -17857.1738
94 -17057.9238 -18973.7988
95 -13329.4238 -17057.9238
96 -15020.9301 -13329.4238
97 -26864.6538 -15020.9301
98 -31224.0982 -26864.6538
99 -27244.2093 -31224.0982
100 -38576.2093 -27244.2093
101 -22673.5967 -38576.2093
102 -25422.4856 -22673.5967
103 -33860.5967 -25422.4856
> 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/7e4w61227791115.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/8vfcs1227791115.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/94gtn1227791115.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/10ozt31227791115.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/11mzov1227791115.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/126u6m1227791115.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/1311yg1227791115.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/14xhov1227791115.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/15k2ht1227791115.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/16rpbc1227791115.tab")
+ }
>
> system("convert tmp/18os01227791115.ps tmp/18os01227791115.png")
> system("convert tmp/229e21227791115.ps tmp/229e21227791115.png")
> system("convert tmp/3iima1227791115.ps tmp/3iima1227791115.png")
> system("convert tmp/4kcny1227791115.ps tmp/4kcny1227791115.png")
> system("convert tmp/5fme81227791115.ps tmp/5fme81227791115.png")
> system("convert tmp/6g4b31227791115.ps tmp/6g4b31227791115.png")
> system("convert tmp/7e4w61227791115.ps tmp/7e4w61227791115.png")
> system("convert tmp/8vfcs1227791115.ps tmp/8vfcs1227791115.png")
> system("convert tmp/94gtn1227791115.ps tmp/94gtn1227791115.png")
> system("convert tmp/10ozt31227791115.ps tmp/10ozt31227791115.png")
>
>
> proc.time()
user system elapsed
6.775 2.853 7.238