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(485,0,464,0,460,0,467,0,460,0,448,0,443,0,436,0,431,0,484,0,510,0,513,0,503,0,471,0,471,0,476,0,475,0,470,0,461,0,455,0,456,0,517,0,525,0,523,0,519,0,509,0,512,0,519,0,517,0,510,0,509,0,501,0,507,0,569,0,580,0,578,0,565,1,547,1,555,1,562,1,561,1,555,1,544,1,537,1,543,1,594,1,611,1,613,1,611,1,594,1,595,1,591,1,589,1,584,1,573,1,567,1,569,1,621,1,629,1,628,1,612,1,595,1,597,1,593,1,590,1,580,1,574,1,573,1,573,1,620,1,626,1,620,1,588,1,566,1,557,1,561,1,549,1,532,1,526,1,511,1,499,1,555,1,565,1,542,1,527,1,510,1,514,1,517,1,508,1,493,1,490,1,469,1,478,1,528,1),dim=c(2,94),dimnames=list(c('Y','D'),1:94))
> y <- array(NA,dim=c(2,94),dimnames=list(c('Y','D'),1:94))
> 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 D M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 485 0 1 0 0 0 0 0 0 0 0 0 0 1
2 464 0 0 1 0 0 0 0 0 0 0 0 0 2
3 460 0 0 0 1 0 0 0 0 0 0 0 0 3
4 467 0 0 0 0 1 0 0 0 0 0 0 0 4
5 460 0 0 0 0 0 1 0 0 0 0 0 0 5
6 448 0 0 0 0 0 0 1 0 0 0 0 0 6
7 443 0 0 0 0 0 0 0 1 0 0 0 0 7
8 436 0 0 0 0 0 0 0 0 1 0 0 0 8
9 431 0 0 0 0 0 0 0 0 0 1 0 0 9
10 484 0 0 0 0 0 0 0 0 0 0 1 0 10
11 510 0 0 0 0 0 0 0 0 0 0 0 1 11
12 513 0 0 0 0 0 0 0 0 0 0 0 0 12
13 503 0 1 0 0 0 0 0 0 0 0 0 0 13
14 471 0 0 1 0 0 0 0 0 0 0 0 0 14
15 471 0 0 0 1 0 0 0 0 0 0 0 0 15
16 476 0 0 0 0 1 0 0 0 0 0 0 0 16
17 475 0 0 0 0 0 1 0 0 0 0 0 0 17
18 470 0 0 0 0 0 0 1 0 0 0 0 0 18
19 461 0 0 0 0 0 0 0 1 0 0 0 0 19
20 455 0 0 0 0 0 0 0 0 1 0 0 0 20
21 456 0 0 0 0 0 0 0 0 0 1 0 0 21
22 517 0 0 0 0 0 0 0 0 0 0 1 0 22
23 525 0 0 0 0 0 0 0 0 0 0 0 1 23
24 523 0 0 0 0 0 0 0 0 0 0 0 0 24
25 519 0 1 0 0 0 0 0 0 0 0 0 0 25
26 509 0 0 1 0 0 0 0 0 0 0 0 0 26
27 512 0 0 0 1 0 0 0 0 0 0 0 0 27
28 519 0 0 0 0 1 0 0 0 0 0 0 0 28
29 517 0 0 0 0 0 1 0 0 0 0 0 0 29
30 510 0 0 0 0 0 0 1 0 0 0 0 0 30
31 509 0 0 0 0 0 0 0 1 0 0 0 0 31
32 501 0 0 0 0 0 0 0 0 1 0 0 0 32
33 507 0 0 0 0 0 0 0 0 0 1 0 0 33
34 569 0 0 0 0 0 0 0 0 0 0 1 0 34
35 580 0 0 0 0 0 0 0 0 0 0 0 1 35
36 578 0 0 0 0 0 0 0 0 0 0 0 0 36
37 565 1 1 0 0 0 0 0 0 0 0 0 0 37
38 547 1 0 1 0 0 0 0 0 0 0 0 0 38
39 555 1 0 0 1 0 0 0 0 0 0 0 0 39
40 562 1 0 0 0 1 0 0 0 0 0 0 0 40
41 561 1 0 0 0 0 1 0 0 0 0 0 0 41
42 555 1 0 0 0 0 0 1 0 0 0 0 0 42
43 544 1 0 0 0 0 0 0 1 0 0 0 0 43
44 537 1 0 0 0 0 0 0 0 1 0 0 0 44
45 543 1 0 0 0 0 0 0 0 0 1 0 0 45
46 594 1 0 0 0 0 0 0 0 0 0 1 0 46
47 611 1 0 0 0 0 0 0 0 0 0 0 1 47
48 613 1 0 0 0 0 0 0 0 0 0 0 0 48
49 611 1 1 0 0 0 0 0 0 0 0 0 0 49
50 594 1 0 1 0 0 0 0 0 0 0 0 0 50
51 595 1 0 0 1 0 0 0 0 0 0 0 0 51
52 591 1 0 0 0 1 0 0 0 0 0 0 0 52
53 589 1 0 0 0 0 1 0 0 0 0 0 0 53
54 584 1 0 0 0 0 0 1 0 0 0 0 0 54
55 573 1 0 0 0 0 0 0 1 0 0 0 0 55
56 567 1 0 0 0 0 0 0 0 1 0 0 0 56
57 569 1 0 0 0 0 0 0 0 0 1 0 0 57
58 621 1 0 0 0 0 0 0 0 0 0 1 0 58
59 629 1 0 0 0 0 0 0 0 0 0 0 1 59
60 628 1 0 0 0 0 0 0 0 0 0 0 0 60
61 612 1 1 0 0 0 0 0 0 0 0 0 0 61
62 595 1 0 1 0 0 0 0 0 0 0 0 0 62
63 597 1 0 0 1 0 0 0 0 0 0 0 0 63
64 593 1 0 0 0 1 0 0 0 0 0 0 0 64
65 590 1 0 0 0 0 1 0 0 0 0 0 0 65
66 580 1 0 0 0 0 0 1 0 0 0 0 0 66
67 574 1 0 0 0 0 0 0 1 0 0 0 0 67
68 573 1 0 0 0 0 0 0 0 1 0 0 0 68
69 573 1 0 0 0 0 0 0 0 0 1 0 0 69
70 620 1 0 0 0 0 0 0 0 0 0 1 0 70
71 626 1 0 0 0 0 0 0 0 0 0 0 1 71
72 620 1 0 0 0 0 0 0 0 0 0 0 0 72
73 588 1 1 0 0 0 0 0 0 0 0 0 0 73
74 566 1 0 1 0 0 0 0 0 0 0 0 0 74
75 557 1 0 0 1 0 0 0 0 0 0 0 0 75
76 561 1 0 0 0 1 0 0 0 0 0 0 0 76
77 549 1 0 0 0 0 1 0 0 0 0 0 0 77
78 532 1 0 0 0 0 0 1 0 0 0 0 0 78
79 526 1 0 0 0 0 0 0 1 0 0 0 0 79
80 511 1 0 0 0 0 0 0 0 1 0 0 0 80
81 499 1 0 0 0 0 0 0 0 0 1 0 0 81
82 555 1 0 0 0 0 0 0 0 0 0 1 0 82
83 565 1 0 0 0 0 0 0 0 0 0 0 1 83
84 542 1 0 0 0 0 0 0 0 0 0 0 0 84
85 527 1 1 0 0 0 0 0 0 0 0 0 0 85
86 510 1 0 1 0 0 0 0 0 0 0 0 0 86
87 514 1 0 0 1 0 0 0 0 0 0 0 0 87
88 517 1 0 0 0 1 0 0 0 0 0 0 0 88
89 508 1 0 0 0 0 1 0 0 0 0 0 0 89
90 493 1 0 0 0 0 0 1 0 0 0 0 0 90
91 490 1 0 0 0 0 0 0 1 0 0 0 0 91
92 469 1 0 0 0 0 0 0 0 1 0 0 0 92
93 478 1 0 0 0 0 0 0 0 0 1 0 0 93
94 528 1 0 0 0 0 0 0 0 0 0 1 0 94
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) D M1 M2 M3 M4
544.9825 99.8638 -30.8935 -49.5562 -48.3439 -44.6316
M5 M6 M7 M8 M9 M10
-48.6693 -57.7070 -63.6197 -71.9074 -70.4451 -15.8578
M11 t
3.5556 -0.5873
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-53.513 -27.678 -6.517 28.075 59.843
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 544.9825 13.5699 40.161 < 2e-16 ***
D 99.8638 12.7472 7.834 1.71e-11 ***
M1 -30.8935 16.6229 -1.858 0.066778 .
M2 -49.5562 16.6007 -2.985 0.003759 **
M3 -48.3439 16.5816 -2.916 0.004605 **
M4 -44.6316 16.5657 -2.694 0.008595 **
M5 -48.6693 16.5529 -2.940 0.004287 **
M6 -57.7070 16.5433 -3.488 0.000793 ***
M7 -63.6197 16.5369 -3.847 0.000239 ***
M8 -71.9074 16.5336 -4.349 3.99e-05 ***
M9 -70.4451 16.5336 -4.261 5.52e-05 ***
M10 -15.8578 16.5367 -0.959 0.340475
M11 3.5556 17.0728 0.208 0.835557
t -0.5873 0.2293 -2.561 0.012306 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 31.94 on 80 degrees of freedom
Multiple R-squared: 0.6771, Adjusted R-squared: 0.6246
F-statistic: 12.9 on 13 and 80 DF, p-value: 9.964e-15
> 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,] 2.870925e-03 5.741850e-03 0.99712908
[2,] 1.273474e-03 2.546948e-03 0.99872653
[3,] 2.336497e-04 4.672995e-04 0.99976635
[4,] 4.559013e-05 9.118026e-05 0.99995441
[5,] 2.471297e-05 4.942593e-05 0.99997529
[6,] 4.992313e-05 9.984625e-05 0.99995008
[7,] 1.297373e-05 2.594747e-05 0.99998703
[8,] 4.901122e-06 9.802245e-06 0.99999510
[9,] 1.023565e-06 2.047129e-06 0.99999898
[10,] 1.920625e-06 3.841249e-06 0.99999808
[11,] 3.991880e-06 7.983760e-06 0.99999601
[12,] 4.518617e-06 9.037233e-06 0.99999548
[13,] 4.466448e-06 8.932896e-06 0.99999553
[14,] 3.893245e-06 7.786489e-06 0.99999611
[15,] 5.455268e-06 1.091054e-05 0.99999454
[16,] 4.573224e-06 9.146449e-06 0.99999543
[17,] 7.551272e-06 1.510254e-05 0.99999245
[18,] 1.487572e-05 2.975144e-05 0.99998512
[19,] 1.190863e-05 2.381726e-05 0.99998809
[20,] 7.191534e-06 1.438307e-05 0.99999281
[21,] 5.576256e-06 1.115251e-05 0.99999442
[22,] 5.230548e-06 1.046110e-05 0.99999477
[23,] 5.489364e-06 1.097873e-05 0.99999451
[24,] 5.260541e-06 1.052108e-05 0.99999474
[25,] 5.228530e-06 1.045706e-05 0.99999477
[26,] 5.285337e-06 1.057067e-05 0.99999471
[27,] 6.975044e-06 1.395009e-05 0.99999302
[28,] 1.096806e-05 2.193612e-05 0.99998903
[29,] 2.021613e-05 4.043226e-05 0.99997978
[30,] 5.827095e-05 1.165419e-04 0.99994173
[31,] 1.796227e-04 3.592454e-04 0.99982038
[32,] 3.880321e-04 7.760641e-04 0.99961197
[33,] 4.269306e-04 8.538612e-04 0.99957307
[34,] 5.504500e-04 1.100900e-03 0.99944955
[35,] 6.851644e-04 1.370329e-03 0.99931484
[36,] 1.165979e-03 2.331959e-03 0.99883402
[37,] 1.733106e-03 3.466211e-03 0.99826689
[38,] 1.838506e-03 3.677012e-03 0.99816149
[39,] 3.725059e-03 7.450117e-03 0.99627494
[40,] 6.362551e-03 1.272510e-02 0.99363745
[41,] 1.187774e-02 2.375549e-02 0.98812226
[42,] 3.208406e-02 6.416812e-02 0.96791594
[43,] 1.074422e-01 2.148843e-01 0.89255783
[44,] 1.605026e-01 3.210052e-01 0.83949739
[45,] 2.308682e-01 4.617364e-01 0.76913182
[46,] 2.634097e-01 5.268194e-01 0.73659032
[47,] 2.619808e-01 5.239615e-01 0.73801923
[48,] 3.703383e-01 7.406766e-01 0.62966168
[49,] 4.093390e-01 8.186780e-01 0.59066099
[50,] 4.182335e-01 8.364670e-01 0.58176650
[51,] 4.132019e-01 8.264038e-01 0.58679811
[52,] 3.833520e-01 7.667039e-01 0.61664803
[53,] 3.823454e-01 7.646908e-01 0.61765459
[54,] 3.760257e-01 7.520514e-01 0.62397429
[55,] 4.351982e-01 8.703965e-01 0.56480175
[56,] 8.365539e-01 3.268922e-01 0.16344612
[57,] 9.457292e-01 1.085416e-01 0.05427079
[58,] 9.803103e-01 3.937932e-02 0.01968966
[59,] 9.736717e-01 5.265650e-02 0.02632825
[60,] 9.616352e-01 7.672952e-02 0.03836476
[61,] 9.287428e-01 1.425144e-01 0.07125718
> postscript(file="/var/www/html/rcomp/tmp/12evk1229938421.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/2jaf81229938421.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/3oha41229938421.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/4ikcj1229938421.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/51lph1229938421.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 = 94
Frequency = 1
1 2 3 4 5 6
-28.5016734 -30.2516734 -34.8766734 -31.0016734 -33.3766734 -35.7516734
7 8 9 10 11 12
-34.2516734 -32.3766734 -38.2516734 -39.2516734 -32.0777276 -24.9348704
13 14 15 16 17 18
-3.4540902 -16.2040902 -16.8290902 -14.9540902 -11.3290902 -6.7040902
19 20 21 22 23 24
-9.2040902 -6.3290902 -6.2040902 0.7959098 -10.0301444 -7.8872872
25 26 27 28 29 30
19.5934930 28.8434930 31.2184930 35.0934930 37.7184930 40.3434930
31 32 33 34 35 36
45.8434930 46.7184930 51.8434930 59.8434930 52.0174388 54.1602959
37 38 39 40 41 42
-27.2227123 -25.9727123 -18.5977123 -14.7227123 -11.0977123 -7.4727123
43 44 45 46 47 48
-11.9727123 -10.0977123 -4.9727123 -7.9727123 -9.7987665 -3.6559094
49 50 51 52 53 54
25.8248709 28.0748709 28.4498709 21.3248709 23.9498709 28.5748709
55 56 57 58 59 60
24.0748709 26.9498709 28.0748709 26.0748709 15.2488167 18.3916738
61 62 63 64 65 66
33.8724541 36.1224541 37.4974541 30.3724541 31.9974541 31.6224541
67 68 69 70 71 72
32.1224541 39.9974541 39.1224541 32.1224541 19.2963999 17.4392570
73 74 75 76 77 78
16.9200373 14.1700373 4.5450373 5.4200373 -1.9549627 -9.3299627
79 80 81 82 83 84
-8.8299627 -14.9549627 -27.8299627 -25.8299627 -34.6560169 -53.5131598
85 86 87 88 89 90
-37.0323795 -34.7823795 -31.4073795 -31.5323795 -35.9073795 -41.2823795
91 92 93 94
-37.7823795 -49.9073795 -41.7823795 -45.7823795
> postscript(file="/var/www/html/rcomp/tmp/6rr661229938421.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 = 94
Frequency = 1
lag(myerror, k = 1) myerror
0 -28.5016734 NA
1 -30.2516734 -28.5016734
2 -34.8766734 -30.2516734
3 -31.0016734 -34.8766734
4 -33.3766734 -31.0016734
5 -35.7516734 -33.3766734
6 -34.2516734 -35.7516734
7 -32.3766734 -34.2516734
8 -38.2516734 -32.3766734
9 -39.2516734 -38.2516734
10 -32.0777276 -39.2516734
11 -24.9348704 -32.0777276
12 -3.4540902 -24.9348704
13 -16.2040902 -3.4540902
14 -16.8290902 -16.2040902
15 -14.9540902 -16.8290902
16 -11.3290902 -14.9540902
17 -6.7040902 -11.3290902
18 -9.2040902 -6.7040902
19 -6.3290902 -9.2040902
20 -6.2040902 -6.3290902
21 0.7959098 -6.2040902
22 -10.0301444 0.7959098
23 -7.8872872 -10.0301444
24 19.5934930 -7.8872872
25 28.8434930 19.5934930
26 31.2184930 28.8434930
27 35.0934930 31.2184930
28 37.7184930 35.0934930
29 40.3434930 37.7184930
30 45.8434930 40.3434930
31 46.7184930 45.8434930
32 51.8434930 46.7184930
33 59.8434930 51.8434930
34 52.0174388 59.8434930
35 54.1602959 52.0174388
36 -27.2227123 54.1602959
37 -25.9727123 -27.2227123
38 -18.5977123 -25.9727123
39 -14.7227123 -18.5977123
40 -11.0977123 -14.7227123
41 -7.4727123 -11.0977123
42 -11.9727123 -7.4727123
43 -10.0977123 -11.9727123
44 -4.9727123 -10.0977123
45 -7.9727123 -4.9727123
46 -9.7987665 -7.9727123
47 -3.6559094 -9.7987665
48 25.8248709 -3.6559094
49 28.0748709 25.8248709
50 28.4498709 28.0748709
51 21.3248709 28.4498709
52 23.9498709 21.3248709
53 28.5748709 23.9498709
54 24.0748709 28.5748709
55 26.9498709 24.0748709
56 28.0748709 26.9498709
57 26.0748709 28.0748709
58 15.2488167 26.0748709
59 18.3916738 15.2488167
60 33.8724541 18.3916738
61 36.1224541 33.8724541
62 37.4974541 36.1224541
63 30.3724541 37.4974541
64 31.9974541 30.3724541
65 31.6224541 31.9974541
66 32.1224541 31.6224541
67 39.9974541 32.1224541
68 39.1224541 39.9974541
69 32.1224541 39.1224541
70 19.2963999 32.1224541
71 17.4392570 19.2963999
72 16.9200373 17.4392570
73 14.1700373 16.9200373
74 4.5450373 14.1700373
75 5.4200373 4.5450373
76 -1.9549627 5.4200373
77 -9.3299627 -1.9549627
78 -8.8299627 -9.3299627
79 -14.9549627 -8.8299627
80 -27.8299627 -14.9549627
81 -25.8299627 -27.8299627
82 -34.6560169 -25.8299627
83 -53.5131598 -34.6560169
84 -37.0323795 -53.5131598
85 -34.7823795 -37.0323795
86 -31.4073795 -34.7823795
87 -31.5323795 -31.4073795
88 -35.9073795 -31.5323795
89 -41.2823795 -35.9073795
90 -37.7823795 -41.2823795
91 -49.9073795 -37.7823795
92 -41.7823795 -49.9073795
93 -45.7823795 -41.7823795
94 NA -45.7823795
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -30.2516734 -28.5016734
[2,] -34.8766734 -30.2516734
[3,] -31.0016734 -34.8766734
[4,] -33.3766734 -31.0016734
[5,] -35.7516734 -33.3766734
[6,] -34.2516734 -35.7516734
[7,] -32.3766734 -34.2516734
[8,] -38.2516734 -32.3766734
[9,] -39.2516734 -38.2516734
[10,] -32.0777276 -39.2516734
[11,] -24.9348704 -32.0777276
[12,] -3.4540902 -24.9348704
[13,] -16.2040902 -3.4540902
[14,] -16.8290902 -16.2040902
[15,] -14.9540902 -16.8290902
[16,] -11.3290902 -14.9540902
[17,] -6.7040902 -11.3290902
[18,] -9.2040902 -6.7040902
[19,] -6.3290902 -9.2040902
[20,] -6.2040902 -6.3290902
[21,] 0.7959098 -6.2040902
[22,] -10.0301444 0.7959098
[23,] -7.8872872 -10.0301444
[24,] 19.5934930 -7.8872872
[25,] 28.8434930 19.5934930
[26,] 31.2184930 28.8434930
[27,] 35.0934930 31.2184930
[28,] 37.7184930 35.0934930
[29,] 40.3434930 37.7184930
[30,] 45.8434930 40.3434930
[31,] 46.7184930 45.8434930
[32,] 51.8434930 46.7184930
[33,] 59.8434930 51.8434930
[34,] 52.0174388 59.8434930
[35,] 54.1602959 52.0174388
[36,] -27.2227123 54.1602959
[37,] -25.9727123 -27.2227123
[38,] -18.5977123 -25.9727123
[39,] -14.7227123 -18.5977123
[40,] -11.0977123 -14.7227123
[41,] -7.4727123 -11.0977123
[42,] -11.9727123 -7.4727123
[43,] -10.0977123 -11.9727123
[44,] -4.9727123 -10.0977123
[45,] -7.9727123 -4.9727123
[46,] -9.7987665 -7.9727123
[47,] -3.6559094 -9.7987665
[48,] 25.8248709 -3.6559094
[49,] 28.0748709 25.8248709
[50,] 28.4498709 28.0748709
[51,] 21.3248709 28.4498709
[52,] 23.9498709 21.3248709
[53,] 28.5748709 23.9498709
[54,] 24.0748709 28.5748709
[55,] 26.9498709 24.0748709
[56,] 28.0748709 26.9498709
[57,] 26.0748709 28.0748709
[58,] 15.2488167 26.0748709
[59,] 18.3916738 15.2488167
[60,] 33.8724541 18.3916738
[61,] 36.1224541 33.8724541
[62,] 37.4974541 36.1224541
[63,] 30.3724541 37.4974541
[64,] 31.9974541 30.3724541
[65,] 31.6224541 31.9974541
[66,] 32.1224541 31.6224541
[67,] 39.9974541 32.1224541
[68,] 39.1224541 39.9974541
[69,] 32.1224541 39.1224541
[70,] 19.2963999 32.1224541
[71,] 17.4392570 19.2963999
[72,] 16.9200373 17.4392570
[73,] 14.1700373 16.9200373
[74,] 4.5450373 14.1700373
[75,] 5.4200373 4.5450373
[76,] -1.9549627 5.4200373
[77,] -9.3299627 -1.9549627
[78,] -8.8299627 -9.3299627
[79,] -14.9549627 -8.8299627
[80,] -27.8299627 -14.9549627
[81,] -25.8299627 -27.8299627
[82,] -34.6560169 -25.8299627
[83,] -53.5131598 -34.6560169
[84,] -37.0323795 -53.5131598
[85,] -34.7823795 -37.0323795
[86,] -31.4073795 -34.7823795
[87,] -31.5323795 -31.4073795
[88,] -35.9073795 -31.5323795
[89,] -41.2823795 -35.9073795
[90,] -37.7823795 -41.2823795
[91,] -49.9073795 -37.7823795
[92,] -41.7823795 -49.9073795
[93,] -45.7823795 -41.7823795
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -30.2516734 -28.5016734
2 -34.8766734 -30.2516734
3 -31.0016734 -34.8766734
4 -33.3766734 -31.0016734
5 -35.7516734 -33.3766734
6 -34.2516734 -35.7516734
7 -32.3766734 -34.2516734
8 -38.2516734 -32.3766734
9 -39.2516734 -38.2516734
10 -32.0777276 -39.2516734
11 -24.9348704 -32.0777276
12 -3.4540902 -24.9348704
13 -16.2040902 -3.4540902
14 -16.8290902 -16.2040902
15 -14.9540902 -16.8290902
16 -11.3290902 -14.9540902
17 -6.7040902 -11.3290902
18 -9.2040902 -6.7040902
19 -6.3290902 -9.2040902
20 -6.2040902 -6.3290902
21 0.7959098 -6.2040902
22 -10.0301444 0.7959098
23 -7.8872872 -10.0301444
24 19.5934930 -7.8872872
25 28.8434930 19.5934930
26 31.2184930 28.8434930
27 35.0934930 31.2184930
28 37.7184930 35.0934930
29 40.3434930 37.7184930
30 45.8434930 40.3434930
31 46.7184930 45.8434930
32 51.8434930 46.7184930
33 59.8434930 51.8434930
34 52.0174388 59.8434930
35 54.1602959 52.0174388
36 -27.2227123 54.1602959
37 -25.9727123 -27.2227123
38 -18.5977123 -25.9727123
39 -14.7227123 -18.5977123
40 -11.0977123 -14.7227123
41 -7.4727123 -11.0977123
42 -11.9727123 -7.4727123
43 -10.0977123 -11.9727123
44 -4.9727123 -10.0977123
45 -7.9727123 -4.9727123
46 -9.7987665 -7.9727123
47 -3.6559094 -9.7987665
48 25.8248709 -3.6559094
49 28.0748709 25.8248709
50 28.4498709 28.0748709
51 21.3248709 28.4498709
52 23.9498709 21.3248709
53 28.5748709 23.9498709
54 24.0748709 28.5748709
55 26.9498709 24.0748709
56 28.0748709 26.9498709
57 26.0748709 28.0748709
58 15.2488167 26.0748709
59 18.3916738 15.2488167
60 33.8724541 18.3916738
61 36.1224541 33.8724541
62 37.4974541 36.1224541
63 30.3724541 37.4974541
64 31.9974541 30.3724541
65 31.6224541 31.9974541
66 32.1224541 31.6224541
67 39.9974541 32.1224541
68 39.1224541 39.9974541
69 32.1224541 39.1224541
70 19.2963999 32.1224541
71 17.4392570 19.2963999
72 16.9200373 17.4392570
73 14.1700373 16.9200373
74 4.5450373 14.1700373
75 5.4200373 4.5450373
76 -1.9549627 5.4200373
77 -9.3299627 -1.9549627
78 -8.8299627 -9.3299627
79 -14.9549627 -8.8299627
80 -27.8299627 -14.9549627
81 -25.8299627 -27.8299627
82 -34.6560169 -25.8299627
83 -53.5131598 -34.6560169
84 -37.0323795 -53.5131598
85 -34.7823795 -37.0323795
86 -31.4073795 -34.7823795
87 -31.5323795 -31.4073795
88 -35.9073795 -31.5323795
89 -41.2823795 -35.9073795
90 -37.7823795 -41.2823795
91 -49.9073795 -37.7823795
92 -41.7823795 -49.9073795
93 -45.7823795 -41.7823795
> 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/7vmco1229938421.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/8hf2q1229938421.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/95akz1229938421.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/10878c1229938421.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/11sk7v1229938421.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/12lssi1229938421.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/13f6191229938421.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/14tw5l1229938421.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/15ssum1229938421.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/1673vi1229938421.tab")
+ }
>
> system("convert tmp/12evk1229938421.ps tmp/12evk1229938421.png")
> system("convert tmp/2jaf81229938421.ps tmp/2jaf81229938421.png")
> system("convert tmp/3oha41229938421.ps tmp/3oha41229938421.png")
> system("convert tmp/4ikcj1229938421.ps tmp/4ikcj1229938421.png")
> system("convert tmp/51lph1229938421.ps tmp/51lph1229938421.png")
> system("convert tmp/6rr661229938421.ps tmp/6rr661229938421.png")
> system("convert tmp/7vmco1229938421.ps tmp/7vmco1229938421.png")
> system("convert tmp/8hf2q1229938421.ps tmp/8hf2q1229938421.png")
> system("convert tmp/95akz1229938421.ps tmp/95akz1229938421.png")
> system("convert tmp/10878c1229938421.ps tmp/10878c1229938421.png")
>
>
> proc.time()
user system elapsed
5.735 2.788 6.115