R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(216234
+ ,562325
+ ,213587
+ ,560854
+ ,209465
+ ,555332
+ ,204045
+ ,543599
+ ,200237
+ ,536662
+ ,203666
+ ,542722
+ ,241476
+ ,593530
+ ,260307
+ ,610763
+ ,243324
+ ,612613
+ ,244460
+ ,611324
+ ,233575
+ ,594167
+ ,237217
+ ,595454
+ ,235243
+ ,590865
+ ,230354
+ ,589379
+ ,227184
+ ,584428
+ ,221678
+ ,573100
+ ,217142
+ ,567456
+ ,219452
+ ,569028
+ ,256446
+ ,620735
+ ,265845
+ ,628884
+ ,248624
+ ,628232
+ ,241114
+ ,612117
+ ,229245
+ ,595404
+ ,231805
+ ,597141
+ ,219277
+ ,593408
+ ,219313
+ ,590072
+ ,212610
+ ,579799
+ ,214771
+ ,574205
+ ,211142
+ ,572775
+ ,211457
+ ,572942
+ ,240048
+ ,619567
+ ,240636
+ ,625809
+ ,230580
+ ,619916
+ ,208795
+ ,587625
+ ,197922
+ ,565742
+ ,194596
+ ,557274
+ ,194581
+ ,560576
+ ,185686
+ ,548854
+ ,178106
+ ,531673
+ ,172608
+ ,525919
+ ,167302
+ ,511038
+ ,168053
+ ,498662
+ ,202300
+ ,555362
+ ,202388
+ ,564591
+ ,182516
+ ,541657
+ ,173476
+ ,527070
+ ,166444
+ ,509846
+ ,171297
+ ,514258
+ ,169701
+ ,516922
+ ,164182
+ ,507561
+ ,161914
+ ,492622
+ ,159612
+ ,490243
+ ,151001
+ ,469357
+ ,158114
+ ,477580
+ ,186530
+ ,528379
+ ,187069
+ ,533590
+ ,174330
+ ,517945
+ ,169362
+ ,506174
+ ,166827
+ ,501866
+ ,178037
+ ,516141
+ ,186412
+ ,528222
+ ,189226
+ ,532638
+ ,191563
+ ,536322
+ ,188906
+ ,536535
+ ,186005
+ ,523597
+ ,195309
+ ,536214
+ ,223532
+ ,586570
+ ,226899
+ ,596594
+ ,214126
+ ,580523)
+ ,dim=c(2
+ ,69)
+ ,dimnames=list(c('Y'
+ ,'X')
+ ,1:69))
> y <- array(NA,dim=c(2,69),dimnames=list(c('Y','X'),1:69))
> 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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 216234 562325 1 0 0 0 0 0 0 0 0 0 0 1
2 213587 560854 0 1 0 0 0 0 0 0 0 0 0 2
3 209465 555332 0 0 1 0 0 0 0 0 0 0 0 3
4 204045 543599 0 0 0 1 0 0 0 0 0 0 0 4
5 200237 536662 0 0 0 0 1 0 0 0 0 0 0 5
6 203666 542722 0 0 0 0 0 1 0 0 0 0 0 6
7 241476 593530 0 0 0 0 0 0 1 0 0 0 0 7
8 260307 610763 0 0 0 0 0 0 0 1 0 0 0 8
9 243324 612613 0 0 0 0 0 0 0 0 1 0 0 9
10 244460 611324 0 0 0 0 0 0 0 0 0 1 0 10
11 233575 594167 0 0 0 0 0 0 0 0 0 0 1 11
12 237217 595454 0 0 0 0 0 0 0 0 0 0 0 12
13 235243 590865 1 0 0 0 0 0 0 0 0 0 0 13
14 230354 589379 0 1 0 0 0 0 0 0 0 0 0 14
15 227184 584428 0 0 1 0 0 0 0 0 0 0 0 15
16 221678 573100 0 0 0 1 0 0 0 0 0 0 0 16
17 217142 567456 0 0 0 0 1 0 0 0 0 0 0 17
18 219452 569028 0 0 0 0 0 1 0 0 0 0 0 18
19 256446 620735 0 0 0 0 0 0 1 0 0 0 0 19
20 265845 628884 0 0 0 0 0 0 0 1 0 0 0 20
21 248624 628232 0 0 0 0 0 0 0 0 1 0 0 21
22 241114 612117 0 0 0 0 0 0 0 0 0 1 0 22
23 229245 595404 0 0 0 0 0 0 0 0 0 0 1 23
24 231805 597141 0 0 0 0 0 0 0 0 0 0 0 24
25 219277 593408 1 0 0 0 0 0 0 0 0 0 0 25
26 219313 590072 0 1 0 0 0 0 0 0 0 0 0 26
27 212610 579799 0 0 1 0 0 0 0 0 0 0 0 27
28 214771 574205 0 0 0 1 0 0 0 0 0 0 0 28
29 211142 572775 0 0 0 0 1 0 0 0 0 0 0 29
30 211457 572942 0 0 0 0 0 1 0 0 0 0 0 30
31 240048 619567 0 0 0 0 0 0 1 0 0 0 0 31
32 240636 625809 0 0 0 0 0 0 0 1 0 0 0 32
33 230580 619916 0 0 0 0 0 0 0 0 1 0 0 33
34 208795 587625 0 0 0 0 0 0 0 0 0 1 0 34
35 197922 565742 0 0 0 0 0 0 0 0 0 0 1 35
36 194596 557274 0 0 0 0 0 0 0 0 0 0 0 36
37 194581 560576 1 0 0 0 0 0 0 0 0 0 0 37
38 185686 548854 0 1 0 0 0 0 0 0 0 0 0 38
39 178106 531673 0 0 1 0 0 0 0 0 0 0 0 39
40 172608 525919 0 0 0 1 0 0 0 0 0 0 0 40
41 167302 511038 0 0 0 0 1 0 0 0 0 0 0 41
42 168053 498662 0 0 0 0 0 1 0 0 0 0 0 42
43 202300 555362 0 0 0 0 0 0 1 0 0 0 0 43
44 202388 564591 0 0 0 0 0 0 0 1 0 0 0 44
45 182516 541657 0 0 0 0 0 0 0 0 1 0 0 45
46 173476 527070 0 0 0 0 0 0 0 0 0 1 0 46
47 166444 509846 0 0 0 0 0 0 0 0 0 0 1 47
48 171297 514258 0 0 0 0 0 0 0 0 0 0 0 48
49 169701 516922 1 0 0 0 0 0 0 0 0 0 0 49
50 164182 507561 0 1 0 0 0 0 0 0 0 0 0 50
51 161914 492622 0 0 1 0 0 0 0 0 0 0 0 51
52 159612 490243 0 0 0 1 0 0 0 0 0 0 0 52
53 151001 469357 0 0 0 0 1 0 0 0 0 0 0 53
54 158114 477580 0 0 0 0 0 1 0 0 0 0 0 54
55 186530 528379 0 0 0 0 0 0 1 0 0 0 0 55
56 187069 533590 0 0 0 0 0 0 0 1 0 0 0 56
57 174330 517945 0 0 0 0 0 0 0 0 1 0 0 57
58 169362 506174 0 0 0 0 0 0 0 0 0 1 0 58
59 166827 501866 0 0 0 0 0 0 0 0 0 0 1 59
60 178037 516141 0 0 0 0 0 0 0 0 0 0 0 60
61 186412 528222 1 0 0 0 0 0 0 0 0 0 0 61
62 189226 532638 0 1 0 0 0 0 0 0 0 0 0 62
63 191563 536322 0 0 1 0 0 0 0 0 0 0 0 63
64 188906 536535 0 0 0 1 0 0 0 0 0 0 0 64
65 186005 523597 0 0 0 0 1 0 0 0 0 0 0 65
66 195309 536214 0 0 0 0 0 1 0 0 0 0 0 66
67 223532 586570 0 0 0 0 0 0 1 0 0 0 0 67
68 226899 596594 0 0 0 0 0 0 0 1 0 0 0 68
69 214126 580523 0 0 0 0 0 0 0 0 1 0 0 69
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
-1.337e+05 6.189e-01 -1.755e+03 -2.352e+03 -6.453e+02 1.419e+02
M5 M6 M7 M8 M9 M10
2.031e+03 4.441e+03 5.371e+03 5.272e+03 -3.329e+03 -3.512e+03
M11 t
-2.366e+03 -2.179e+02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-11272 -5354 1170 4315 12483
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.337e+05 1.898e+04 -7.047 3.16e-09 ***
X 6.189e-01 3.119e-02 19.845 < 2e-16 ***
M1 -1.755e+03 4.023e+03 -0.436 0.664348
M2 -2.352e+03 4.024e+03 -0.584 0.561275
M3 -6.453e+02 4.039e+03 -0.160 0.873654
M4 1.419e+02 4.056e+03 0.035 0.972223
M5 2.031e+03 4.105e+03 0.495 0.622775
M6 4.441e+03 4.081e+03 1.088 0.281258
M7 5.371e+03 4.119e+03 1.304 0.197643
M8 5.272e+03 4.204e+03 1.254 0.215168
M9 -3.329e+03 4.133e+03 -0.805 0.424029
M10 -3.512e+03 4.209e+03 -0.835 0.407593
M11 -2.366e+03 4.197e+03 -0.564 0.575195
t -2.179e+02 5.497e+01 -3.963 0.000215 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6633 on 55 degrees of freedom
Multiple R-squared: 0.9565, Adjusted R-squared: 0.9462
F-statistic: 92.96 on 13 and 55 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.0016792082 3.358416e-03 9.983208e-01
[2,] 0.0005851654 1.170331e-03 9.994148e-01
[3,] 0.0002273085 4.546170e-04 9.997727e-01
[4,] 0.0003605789 7.211577e-04 9.996394e-01
[5,] 0.0001340748 2.681496e-04 9.998659e-01
[6,] 0.0004719706 9.439412e-04 9.995280e-01
[7,] 0.0002713509 5.427018e-04 9.997286e-01
[8,] 0.0002908141 5.816282e-04 9.997092e-01
[9,] 0.1372080679 2.744161e-01 8.627919e-01
[10,] 0.1402009414 2.804019e-01 8.597991e-01
[11,] 0.1019787576 2.039575e-01 8.980212e-01
[12,] 0.2721563827 5.443128e-01 7.278436e-01
[13,] 0.3060246810 6.120494e-01 6.939753e-01
[14,] 0.2528194619 5.056389e-01 7.471805e-01
[15,] 0.4298951583 8.597903e-01 5.701048e-01
[16,] 0.9349190669 1.301619e-01 6.508093e-02
[17,] 0.9595718565 8.085629e-02 4.042814e-02
[18,] 0.9634693378 7.306132e-02 3.653066e-02
[19,] 0.9804826844 3.903463e-02 1.951732e-02
[20,] 0.9855968230 2.880635e-02 1.440318e-02
[21,] 0.9871243979 2.575120e-02 1.287560e-02
[22,] 0.9838389249 3.232215e-02 1.616108e-02
[23,] 0.9867114435 2.657711e-02 1.328856e-02
[24,] 0.9787075466 4.258491e-02 2.129245e-02
[25,] 0.9811249815 3.775004e-02 1.887502e-02
[26,] 0.9924358051 1.512839e-02 7.564195e-03
[27,] 0.9982427689 3.514462e-03 1.757231e-03
[28,] 0.9992905151 1.418970e-03 7.094849e-04
[29,] 0.9991296049 1.740790e-03 8.703951e-04
[30,] 0.9977022795 4.595441e-03 2.297721e-03
[31,] 0.9983582836 3.283433e-03 1.641716e-03
[32,] 0.9997417507 5.164986e-04 2.582493e-04
[33,] 0.9990272172 1.945566e-03 9.727828e-04
[34,] 0.9999768353 4.632932e-05 2.316466e-05
[35,] 0.9999932303 1.353931e-05 6.769655e-06
[36,] 0.9999766772 4.664558e-05 2.332279e-05
> postscript(file="/var/www/html/rcomp/tmp/15ccb1258572771.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/20yt21258572771.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/3s7wj1258572771.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/4flzk1258572771.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/54s9p1258572771.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 = 69
Frequency = 1
1 2 3 4 5 6
3891.8472 2970.1756 776.8573 2049.5648 864.0950 -1650.3618
7 8 9 10 11 12
4000.9472 12482.9766 3173.8285 5508.5595 4314.5493 5011.4915
13 14 15 16 17 18
7850.8445 4696.4569 3101.7266 4037.7653 1324.0123 468.3367
19 20 21 22 23 24
4747.2229 9419.6596 1421.0872 4286.1513 1833.3339 1169.7552
25 26 27 28 29 30
-7074.7003 -4159.0578 -5992.8150 -938.7507 -5353.6965 -7334.7682
31 32 33 34 35 36
-8313.4534 -11271.7074 -8861.4404 -10259.4808 -8516.4033 -8749.7485
37 38 39 40 41 42
-8835.4130 -9660.3801 -8095.5331 -10601.4391 -8368.0878 -2149.8555
43 44 45 46 47 48
-3708.3126 -9015.3260 -5873.7846 -5484.4535 -2783.9950 -2810.2253
49 50 51 52 53 54
-4082.0091 -2992.2822 2496.9121 1098.0999 3743.1570 3573.9434
55 56 57 58 59 60
-163.1772 -2532.3090 3230.8139 5949.2235 5152.5151 5378.7271
61 62 63 64 65 66
8249.4307 9145.0876 7712.8520 4354.7598 7790.5200 7092.7053
67 68 69
3436.7730 916.7061 6909.4954
> postscript(file="/var/www/html/rcomp/tmp/6jndk1258572771.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 = 69
Frequency = 1
lag(myerror, k = 1) myerror
0 3891.8472 NA
1 2970.1756 3891.8472
2 776.8573 2970.1756
3 2049.5648 776.8573
4 864.0950 2049.5648
5 -1650.3618 864.0950
6 4000.9472 -1650.3618
7 12482.9766 4000.9472
8 3173.8285 12482.9766
9 5508.5595 3173.8285
10 4314.5493 5508.5595
11 5011.4915 4314.5493
12 7850.8445 5011.4915
13 4696.4569 7850.8445
14 3101.7266 4696.4569
15 4037.7653 3101.7266
16 1324.0123 4037.7653
17 468.3367 1324.0123
18 4747.2229 468.3367
19 9419.6596 4747.2229
20 1421.0872 9419.6596
21 4286.1513 1421.0872
22 1833.3339 4286.1513
23 1169.7552 1833.3339
24 -7074.7003 1169.7552
25 -4159.0578 -7074.7003
26 -5992.8150 -4159.0578
27 -938.7507 -5992.8150
28 -5353.6965 -938.7507
29 -7334.7682 -5353.6965
30 -8313.4534 -7334.7682
31 -11271.7074 -8313.4534
32 -8861.4404 -11271.7074
33 -10259.4808 -8861.4404
34 -8516.4033 -10259.4808
35 -8749.7485 -8516.4033
36 -8835.4130 -8749.7485
37 -9660.3801 -8835.4130
38 -8095.5331 -9660.3801
39 -10601.4391 -8095.5331
40 -8368.0878 -10601.4391
41 -2149.8555 -8368.0878
42 -3708.3126 -2149.8555
43 -9015.3260 -3708.3126
44 -5873.7846 -9015.3260
45 -5484.4535 -5873.7846
46 -2783.9950 -5484.4535
47 -2810.2253 -2783.9950
48 -4082.0091 -2810.2253
49 -2992.2822 -4082.0091
50 2496.9121 -2992.2822
51 1098.0999 2496.9121
52 3743.1570 1098.0999
53 3573.9434 3743.1570
54 -163.1772 3573.9434
55 -2532.3090 -163.1772
56 3230.8139 -2532.3090
57 5949.2235 3230.8139
58 5152.5151 5949.2235
59 5378.7271 5152.5151
60 8249.4307 5378.7271
61 9145.0876 8249.4307
62 7712.8520 9145.0876
63 4354.7598 7712.8520
64 7790.5200 4354.7598
65 7092.7053 7790.5200
66 3436.7730 7092.7053
67 916.7061 3436.7730
68 6909.4954 916.7061
69 NA 6909.4954
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 2970.1756 3891.8472
[2,] 776.8573 2970.1756
[3,] 2049.5648 776.8573
[4,] 864.0950 2049.5648
[5,] -1650.3618 864.0950
[6,] 4000.9472 -1650.3618
[7,] 12482.9766 4000.9472
[8,] 3173.8285 12482.9766
[9,] 5508.5595 3173.8285
[10,] 4314.5493 5508.5595
[11,] 5011.4915 4314.5493
[12,] 7850.8445 5011.4915
[13,] 4696.4569 7850.8445
[14,] 3101.7266 4696.4569
[15,] 4037.7653 3101.7266
[16,] 1324.0123 4037.7653
[17,] 468.3367 1324.0123
[18,] 4747.2229 468.3367
[19,] 9419.6596 4747.2229
[20,] 1421.0872 9419.6596
[21,] 4286.1513 1421.0872
[22,] 1833.3339 4286.1513
[23,] 1169.7552 1833.3339
[24,] -7074.7003 1169.7552
[25,] -4159.0578 -7074.7003
[26,] -5992.8150 -4159.0578
[27,] -938.7507 -5992.8150
[28,] -5353.6965 -938.7507
[29,] -7334.7682 -5353.6965
[30,] -8313.4534 -7334.7682
[31,] -11271.7074 -8313.4534
[32,] -8861.4404 -11271.7074
[33,] -10259.4808 -8861.4404
[34,] -8516.4033 -10259.4808
[35,] -8749.7485 -8516.4033
[36,] -8835.4130 -8749.7485
[37,] -9660.3801 -8835.4130
[38,] -8095.5331 -9660.3801
[39,] -10601.4391 -8095.5331
[40,] -8368.0878 -10601.4391
[41,] -2149.8555 -8368.0878
[42,] -3708.3126 -2149.8555
[43,] -9015.3260 -3708.3126
[44,] -5873.7846 -9015.3260
[45,] -5484.4535 -5873.7846
[46,] -2783.9950 -5484.4535
[47,] -2810.2253 -2783.9950
[48,] -4082.0091 -2810.2253
[49,] -2992.2822 -4082.0091
[50,] 2496.9121 -2992.2822
[51,] 1098.0999 2496.9121
[52,] 3743.1570 1098.0999
[53,] 3573.9434 3743.1570
[54,] -163.1772 3573.9434
[55,] -2532.3090 -163.1772
[56,] 3230.8139 -2532.3090
[57,] 5949.2235 3230.8139
[58,] 5152.5151 5949.2235
[59,] 5378.7271 5152.5151
[60,] 8249.4307 5378.7271
[61,] 9145.0876 8249.4307
[62,] 7712.8520 9145.0876
[63,] 4354.7598 7712.8520
[64,] 7790.5200 4354.7598
[65,] 7092.7053 7790.5200
[66,] 3436.7730 7092.7053
[67,] 916.7061 3436.7730
[68,] 6909.4954 916.7061
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 2970.1756 3891.8472
2 776.8573 2970.1756
3 2049.5648 776.8573
4 864.0950 2049.5648
5 -1650.3618 864.0950
6 4000.9472 -1650.3618
7 12482.9766 4000.9472
8 3173.8285 12482.9766
9 5508.5595 3173.8285
10 4314.5493 5508.5595
11 5011.4915 4314.5493
12 7850.8445 5011.4915
13 4696.4569 7850.8445
14 3101.7266 4696.4569
15 4037.7653 3101.7266
16 1324.0123 4037.7653
17 468.3367 1324.0123
18 4747.2229 468.3367
19 9419.6596 4747.2229
20 1421.0872 9419.6596
21 4286.1513 1421.0872
22 1833.3339 4286.1513
23 1169.7552 1833.3339
24 -7074.7003 1169.7552
25 -4159.0578 -7074.7003
26 -5992.8150 -4159.0578
27 -938.7507 -5992.8150
28 -5353.6965 -938.7507
29 -7334.7682 -5353.6965
30 -8313.4534 -7334.7682
31 -11271.7074 -8313.4534
32 -8861.4404 -11271.7074
33 -10259.4808 -8861.4404
34 -8516.4033 -10259.4808
35 -8749.7485 -8516.4033
36 -8835.4130 -8749.7485
37 -9660.3801 -8835.4130
38 -8095.5331 -9660.3801
39 -10601.4391 -8095.5331
40 -8368.0878 -10601.4391
41 -2149.8555 -8368.0878
42 -3708.3126 -2149.8555
43 -9015.3260 -3708.3126
44 -5873.7846 -9015.3260
45 -5484.4535 -5873.7846
46 -2783.9950 -5484.4535
47 -2810.2253 -2783.9950
48 -4082.0091 -2810.2253
49 -2992.2822 -4082.0091
50 2496.9121 -2992.2822
51 1098.0999 2496.9121
52 3743.1570 1098.0999
53 3573.9434 3743.1570
54 -163.1772 3573.9434
55 -2532.3090 -163.1772
56 3230.8139 -2532.3090
57 5949.2235 3230.8139
58 5152.5151 5949.2235
59 5378.7271 5152.5151
60 8249.4307 5378.7271
61 9145.0876 8249.4307
62 7712.8520 9145.0876
63 4354.7598 7712.8520
64 7790.5200 4354.7598
65 7092.7053 7790.5200
66 3436.7730 7092.7053
67 916.7061 3436.7730
68 6909.4954 916.7061
> 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/7r9p71258572771.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/8erbj1258572771.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/9fx471258572771.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/10hmrm1258572771.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/117osl1258572771.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/126yfq1258572771.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/13baml1258572771.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/14xmap1258572772.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/15k0fi1258572772.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/16hpwa1258572772.tab")
+ }
>
> system("convert tmp/15ccb1258572771.ps tmp/15ccb1258572771.png")
> system("convert tmp/20yt21258572771.ps tmp/20yt21258572771.png")
> system("convert tmp/3s7wj1258572771.ps tmp/3s7wj1258572771.png")
> system("convert tmp/4flzk1258572771.ps tmp/4flzk1258572771.png")
> system("convert tmp/54s9p1258572771.ps tmp/54s9p1258572771.png")
> system("convert tmp/6jndk1258572771.ps tmp/6jndk1258572771.png")
> system("convert tmp/7r9p71258572771.ps tmp/7r9p71258572771.png")
> system("convert tmp/8erbj1258572771.ps tmp/8erbj1258572771.png")
> system("convert tmp/9fx471258572771.ps tmp/9fx471258572771.png")
> system("convert tmp/10hmrm1258572771.ps tmp/10hmrm1258572771.png")
>
>
> proc.time()
user system elapsed
2.584 1.598 4.895