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(519,97.4,517,97,510,105.4,509,102.7,501,98.1,507,104.5,569,87.4,580,89.9,578,109.8,565,111.7,547,98.6,555,96.9,562,95.1,561,97,555,112.7,544,102.9,537,97.4,543,111.4,594,87.4,611,96.8,613,114.1,611,110.3,594,103.9,595,101.6,591,94.6,589,95.9,584,104.7,573,102.8,567,98.1,569,113.9,621,80.9,629,95.7,628,113.2,612,105.9,595,108.8,597,102.3,593,99,590,100.7,580,115.5,574,100.7,573,109.9,573,114.6,620,85.4,626,100.5,620,114.8,588,116.5,566,112.9,557,102,561,106,549,105.3,532,118.8,526,106.1,511,109.3,499,117.2,555,92.5,565,104.2,542,112.5,527,122.4,510,113.3,514,100,517,110.7,508,112.8,493,109.8,490,117.3,469,109.1,478,115.9,528,96,534,99.8,518,116.8,506,115.7,502,99.4,516,94.3),dim=c(2,72),dimnames=list(c('Y','X'),1:72))
> y <- array(NA,dim=c(2,72),dimnames=list(c('Y','X'),1:72))
> 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 519 97.4 1 0 0 0 0 0 0 0 0 0 0 1
2 517 97.0 0 1 0 0 0 0 0 0 0 0 0 2
3 510 105.4 0 0 1 0 0 0 0 0 0 0 0 3
4 509 102.7 0 0 0 1 0 0 0 0 0 0 0 4
5 501 98.1 0 0 0 0 1 0 0 0 0 0 0 5
6 507 104.5 0 0 0 0 0 1 0 0 0 0 0 6
7 569 87.4 0 0 0 0 0 0 1 0 0 0 0 7
8 580 89.9 0 0 0 0 0 0 0 1 0 0 0 8
9 578 109.8 0 0 0 0 0 0 0 0 1 0 0 9
10 565 111.7 0 0 0 0 0 0 0 0 0 1 0 10
11 547 98.6 0 0 0 0 0 0 0 0 0 0 1 11
12 555 96.9 0 0 0 0 0 0 0 0 0 0 0 12
13 562 95.1 1 0 0 0 0 0 0 0 0 0 0 13
14 561 97.0 0 1 0 0 0 0 0 0 0 0 0 14
15 555 112.7 0 0 1 0 0 0 0 0 0 0 0 15
16 544 102.9 0 0 0 1 0 0 0 0 0 0 0 16
17 537 97.4 0 0 0 0 1 0 0 0 0 0 0 17
18 543 111.4 0 0 0 0 0 1 0 0 0 0 0 18
19 594 87.4 0 0 0 0 0 0 1 0 0 0 0 19
20 611 96.8 0 0 0 0 0 0 0 1 0 0 0 20
21 613 114.1 0 0 0 0 0 0 0 0 1 0 0 21
22 611 110.3 0 0 0 0 0 0 0 0 0 1 0 22
23 594 103.9 0 0 0 0 0 0 0 0 0 0 1 23
24 595 101.6 0 0 0 0 0 0 0 0 0 0 0 24
25 591 94.6 1 0 0 0 0 0 0 0 0 0 0 25
26 589 95.9 0 1 0 0 0 0 0 0 0 0 0 26
27 584 104.7 0 0 1 0 0 0 0 0 0 0 0 27
28 573 102.8 0 0 0 1 0 0 0 0 0 0 0 28
29 567 98.1 0 0 0 0 1 0 0 0 0 0 0 29
30 569 113.9 0 0 0 0 0 1 0 0 0 0 0 30
31 621 80.9 0 0 0 0 0 0 1 0 0 0 0 31
32 629 95.7 0 0 0 0 0 0 0 1 0 0 0 32
33 628 113.2 0 0 0 0 0 0 0 0 1 0 0 33
34 612 105.9 0 0 0 0 0 0 0 0 0 1 0 34
35 595 108.8 0 0 0 0 0 0 0 0 0 0 1 35
36 597 102.3 0 0 0 0 0 0 0 0 0 0 0 36
37 593 99.0 1 0 0 0 0 0 0 0 0 0 0 37
38 590 100.7 0 1 0 0 0 0 0 0 0 0 0 38
39 580 115.5 0 0 1 0 0 0 0 0 0 0 0 39
40 574 100.7 0 0 0 1 0 0 0 0 0 0 0 40
41 573 109.9 0 0 0 0 1 0 0 0 0 0 0 41
42 573 114.6 0 0 0 0 0 1 0 0 0 0 0 42
43 620 85.4 0 0 0 0 0 0 1 0 0 0 0 43
44 626 100.5 0 0 0 0 0 0 0 1 0 0 0 44
45 620 114.8 0 0 0 0 0 0 0 0 1 0 0 45
46 588 116.5 0 0 0 0 0 0 0 0 0 1 0 46
47 566 112.9 0 0 0 0 0 0 0 0 0 0 1 47
48 557 102.0 0 0 0 0 0 0 0 0 0 0 0 48
49 561 106.0 1 0 0 0 0 0 0 0 0 0 0 49
50 549 105.3 0 1 0 0 0 0 0 0 0 0 0 50
51 532 118.8 0 0 1 0 0 0 0 0 0 0 0 51
52 526 106.1 0 0 0 1 0 0 0 0 0 0 0 52
53 511 109.3 0 0 0 0 1 0 0 0 0 0 0 53
54 499 117.2 0 0 0 0 0 1 0 0 0 0 0 54
55 555 92.5 0 0 0 0 0 0 1 0 0 0 0 55
56 565 104.2 0 0 0 0 0 0 0 1 0 0 0 56
57 542 112.5 0 0 0 0 0 0 0 0 1 0 0 57
58 527 122.4 0 0 0 0 0 0 0 0 0 1 0 58
59 510 113.3 0 0 0 0 0 0 0 0 0 0 1 59
60 514 100.0 0 0 0 0 0 0 0 0 0 0 0 60
61 517 110.7 1 0 0 0 0 0 0 0 0 0 0 61
62 508 112.8 0 1 0 0 0 0 0 0 0 0 0 62
63 493 109.8 0 0 1 0 0 0 0 0 0 0 0 63
64 490 117.3 0 0 0 1 0 0 0 0 0 0 0 64
65 469 109.1 0 0 0 0 1 0 0 0 0 0 0 65
66 478 115.9 0 0 0 0 0 1 0 0 0 0 0 66
67 528 96.0 0 0 0 0 0 0 1 0 0 0 0 67
68 534 99.8 0 0 0 0 0 0 0 1 0 0 0 68
69 518 116.8 0 0 0 0 0 0 0 0 1 0 0 69
70 506 115.7 0 0 0 0 0 0 0 0 0 1 0 70
71 502 99.4 0 0 0 0 0 0 0 0 0 0 1 71
72 516 94.3 0 0 0 0 0 0 0 0 0 0 0 72
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
611.4093 -0.2702 -5.7988 -9.6796 -16.3713 -23.5672
M5 M6 M7 M8 M9 M10
-33.0244 -27.9999 19.0254 31.9598 29.2274 14.9728
M11 t
-2.2276 -0.6869
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-59.6015 -28.2166 -0.6037 32.4971 52.4769
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 611.4093 102.1328 5.986 1.44e-07 ***
X -0.2702 1.0828 -0.250 0.8038
M1 -5.7988 20.7558 -0.279 0.7809
M2 -9.6796 20.8779 -0.464 0.6447
M3 -16.3713 24.8951 -0.658 0.5134
M4 -23.5672 21.9068 -1.076 0.2865
M5 -33.0244 21.2479 -1.554 0.1256
M6 -27.9999 25.6848 -1.090 0.2802
M7 19.0254 23.3956 0.813 0.4194
M8 31.9598 20.4844 1.560 0.1242
M9 29.2274 25.7640 1.134 0.2613
M10 14.9728 25.8011 0.580 0.5639
M11 -2.2276 21.7172 -0.103 0.9187
t -0.6869 0.2633 -2.608 0.0116 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 35.39 on 58 degrees of freedom
Multiple R-squared: 0.3884, Adjusted R-squared: 0.2513
F-statistic: 2.834 on 13 and 58 DF, p-value: 0.003305
> 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,] 6.447781e-03 1.289556e-02 0.993552219
[2,] 2.252797e-03 4.505593e-03 0.997747203
[3,] 5.442590e-03 1.088518e-02 0.994557410
[4,] 3.223540e-03 6.447081e-03 0.996776460
[5,] 1.234663e-03 2.469326e-03 0.998765337
[6,] 7.435268e-04 1.487054e-03 0.999256473
[7,] 5.988976e-04 1.197795e-03 0.999401102
[8,] 2.741358e-04 5.482716e-04 0.999725864
[9,] 2.041542e-04 4.083084e-04 0.999795846
[10,] 1.291758e-04 2.583517e-04 0.999870824
[11,] 5.472841e-05 1.094568e-04 0.999945272
[12,] 5.942502e-05 1.188500e-04 0.999940575
[13,] 4.821973e-05 9.643947e-05 0.999951780
[14,] 8.122248e-05 1.624450e-04 0.999918778
[15,] 1.548687e-04 3.097373e-04 0.999845131
[16,] 1.470859e-03 2.941719e-03 0.998529141
[17,] 3.464540e-03 6.929080e-03 0.996535460
[18,] 1.264041e-02 2.528082e-02 0.987359589
[19,] 4.121993e-02 8.243987e-02 0.958780066
[20,] 7.157121e-02 1.431424e-01 0.928428788
[21,] 1.216542e-01 2.433084e-01 0.878345778
[22,] 1.504507e-01 3.009014e-01 0.849549306
[23,] 1.679581e-01 3.359162e-01 0.832041896
[24,] 1.894032e-01 3.788065e-01 0.810596768
[25,] 2.076902e-01 4.153805e-01 0.792309752
[26,] 2.536288e-01 5.072577e-01 0.746371173
[27,] 2.744597e-01 5.489195e-01 0.725540268
[28,] 3.533508e-01 7.067016e-01 0.646649216
[29,] 7.827502e-01 4.344995e-01 0.217249766
[30,] 9.325218e-01 1.349563e-01 0.067478164
[31,] 9.815722e-01 3.685558e-02 0.018427792
[32,] 9.891223e-01 2.175530e-02 0.010877651
[33,] 9.878573e-01 2.428531e-02 0.012142657
[34,] 9.844886e-01 3.102281e-02 0.015511406
[35,] 9.931448e-01 1.371038e-02 0.006855189
[36,] 9.872704e-01 2.545926e-02 0.012729632
[37,] 9.906325e-01 1.873505e-02 0.009367527
[38,] 9.771745e-01 4.565105e-02 0.022825527
[39,] 9.368677e-01 1.262647e-01 0.063132348
> postscript(file="/var/www/html/rcomp/tmp/11f9e1258624273.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/2kd3z1258624273.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/3wwu11258624273.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/4idn71258624273.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/5pkv31258624273.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 = 72
Frequency = 1
1 2 3 4 5 6
-59.6015484 -57.1420569 -54.4933782 -48.3402952 -47.4393288 -44.0473706
7 8 9 10 11 12
-33.0070046 -33.5789137 -26.7817131 -24.3267971 -27.9798231 -21.9799371
13 14 15 16 17 18
-8.9806673 -4.8996074 0.7218755 -5.0437963 -3.3860523 2.0597842
19 20 21 22 23 24
0.2354449 7.5282411 17.6227991 29.5373064 28.6949363 27.5326739
25 26 27 28 29 30
28.1266587 31.0455703 35.8023478 32.1716285 35.0455703 36.9778516
31 32 33 34 35 36
33.7212879 33.4734188 40.6220262 37.5906685 39.2615968 37.9642965
37 38 39 40 41 42
39.5581956 41.5852061 42.9634665 40.8465591 52.4769361 49.4094741
43 44 45 46 47 48
42.1798496 40.0130546 41.2968711 24.6977377 19.6120596 6.1256718
49 50 51 52 53 54
17.6923751 10.0707924 4.0977316 2.5483432 -1.4427627 -15.6454338
55 56 57 58 59 60
-12.6589462 -11.7445815 -29.0822478 -26.4653546 -28.0373920 -29.1723729
61 62 63 64 65 66
-16.7950138 -20.6599045 -29.0920432 -22.1824393 -35.2543626 -28.7543056
67 68 69 70 71 72
-30.4706316 -35.6912194 -43.6777355 -41.0335609 -31.5513778 -20.4703321
> postscript(file="/var/www/html/rcomp/tmp/6d1b21258624273.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 = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 -59.6015484 NA
1 -57.1420569 -59.6015484
2 -54.4933782 -57.1420569
3 -48.3402952 -54.4933782
4 -47.4393288 -48.3402952
5 -44.0473706 -47.4393288
6 -33.0070046 -44.0473706
7 -33.5789137 -33.0070046
8 -26.7817131 -33.5789137
9 -24.3267971 -26.7817131
10 -27.9798231 -24.3267971
11 -21.9799371 -27.9798231
12 -8.9806673 -21.9799371
13 -4.8996074 -8.9806673
14 0.7218755 -4.8996074
15 -5.0437963 0.7218755
16 -3.3860523 -5.0437963
17 2.0597842 -3.3860523
18 0.2354449 2.0597842
19 7.5282411 0.2354449
20 17.6227991 7.5282411
21 29.5373064 17.6227991
22 28.6949363 29.5373064
23 27.5326739 28.6949363
24 28.1266587 27.5326739
25 31.0455703 28.1266587
26 35.8023478 31.0455703
27 32.1716285 35.8023478
28 35.0455703 32.1716285
29 36.9778516 35.0455703
30 33.7212879 36.9778516
31 33.4734188 33.7212879
32 40.6220262 33.4734188
33 37.5906685 40.6220262
34 39.2615968 37.5906685
35 37.9642965 39.2615968
36 39.5581956 37.9642965
37 41.5852061 39.5581956
38 42.9634665 41.5852061
39 40.8465591 42.9634665
40 52.4769361 40.8465591
41 49.4094741 52.4769361
42 42.1798496 49.4094741
43 40.0130546 42.1798496
44 41.2968711 40.0130546
45 24.6977377 41.2968711
46 19.6120596 24.6977377
47 6.1256718 19.6120596
48 17.6923751 6.1256718
49 10.0707924 17.6923751
50 4.0977316 10.0707924
51 2.5483432 4.0977316
52 -1.4427627 2.5483432
53 -15.6454338 -1.4427627
54 -12.6589462 -15.6454338
55 -11.7445815 -12.6589462
56 -29.0822478 -11.7445815
57 -26.4653546 -29.0822478
58 -28.0373920 -26.4653546
59 -29.1723729 -28.0373920
60 -16.7950138 -29.1723729
61 -20.6599045 -16.7950138
62 -29.0920432 -20.6599045
63 -22.1824393 -29.0920432
64 -35.2543626 -22.1824393
65 -28.7543056 -35.2543626
66 -30.4706316 -28.7543056
67 -35.6912194 -30.4706316
68 -43.6777355 -35.6912194
69 -41.0335609 -43.6777355
70 -31.5513778 -41.0335609
71 -20.4703321 -31.5513778
72 NA -20.4703321
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -57.1420569 -59.6015484
[2,] -54.4933782 -57.1420569
[3,] -48.3402952 -54.4933782
[4,] -47.4393288 -48.3402952
[5,] -44.0473706 -47.4393288
[6,] -33.0070046 -44.0473706
[7,] -33.5789137 -33.0070046
[8,] -26.7817131 -33.5789137
[9,] -24.3267971 -26.7817131
[10,] -27.9798231 -24.3267971
[11,] -21.9799371 -27.9798231
[12,] -8.9806673 -21.9799371
[13,] -4.8996074 -8.9806673
[14,] 0.7218755 -4.8996074
[15,] -5.0437963 0.7218755
[16,] -3.3860523 -5.0437963
[17,] 2.0597842 -3.3860523
[18,] 0.2354449 2.0597842
[19,] 7.5282411 0.2354449
[20,] 17.6227991 7.5282411
[21,] 29.5373064 17.6227991
[22,] 28.6949363 29.5373064
[23,] 27.5326739 28.6949363
[24,] 28.1266587 27.5326739
[25,] 31.0455703 28.1266587
[26,] 35.8023478 31.0455703
[27,] 32.1716285 35.8023478
[28,] 35.0455703 32.1716285
[29,] 36.9778516 35.0455703
[30,] 33.7212879 36.9778516
[31,] 33.4734188 33.7212879
[32,] 40.6220262 33.4734188
[33,] 37.5906685 40.6220262
[34,] 39.2615968 37.5906685
[35,] 37.9642965 39.2615968
[36,] 39.5581956 37.9642965
[37,] 41.5852061 39.5581956
[38,] 42.9634665 41.5852061
[39,] 40.8465591 42.9634665
[40,] 52.4769361 40.8465591
[41,] 49.4094741 52.4769361
[42,] 42.1798496 49.4094741
[43,] 40.0130546 42.1798496
[44,] 41.2968711 40.0130546
[45,] 24.6977377 41.2968711
[46,] 19.6120596 24.6977377
[47,] 6.1256718 19.6120596
[48,] 17.6923751 6.1256718
[49,] 10.0707924 17.6923751
[50,] 4.0977316 10.0707924
[51,] 2.5483432 4.0977316
[52,] -1.4427627 2.5483432
[53,] -15.6454338 -1.4427627
[54,] -12.6589462 -15.6454338
[55,] -11.7445815 -12.6589462
[56,] -29.0822478 -11.7445815
[57,] -26.4653546 -29.0822478
[58,] -28.0373920 -26.4653546
[59,] -29.1723729 -28.0373920
[60,] -16.7950138 -29.1723729
[61,] -20.6599045 -16.7950138
[62,] -29.0920432 -20.6599045
[63,] -22.1824393 -29.0920432
[64,] -35.2543626 -22.1824393
[65,] -28.7543056 -35.2543626
[66,] -30.4706316 -28.7543056
[67,] -35.6912194 -30.4706316
[68,] -43.6777355 -35.6912194
[69,] -41.0335609 -43.6777355
[70,] -31.5513778 -41.0335609
[71,] -20.4703321 -31.5513778
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -57.1420569 -59.6015484
2 -54.4933782 -57.1420569
3 -48.3402952 -54.4933782
4 -47.4393288 -48.3402952
5 -44.0473706 -47.4393288
6 -33.0070046 -44.0473706
7 -33.5789137 -33.0070046
8 -26.7817131 -33.5789137
9 -24.3267971 -26.7817131
10 -27.9798231 -24.3267971
11 -21.9799371 -27.9798231
12 -8.9806673 -21.9799371
13 -4.8996074 -8.9806673
14 0.7218755 -4.8996074
15 -5.0437963 0.7218755
16 -3.3860523 -5.0437963
17 2.0597842 -3.3860523
18 0.2354449 2.0597842
19 7.5282411 0.2354449
20 17.6227991 7.5282411
21 29.5373064 17.6227991
22 28.6949363 29.5373064
23 27.5326739 28.6949363
24 28.1266587 27.5326739
25 31.0455703 28.1266587
26 35.8023478 31.0455703
27 32.1716285 35.8023478
28 35.0455703 32.1716285
29 36.9778516 35.0455703
30 33.7212879 36.9778516
31 33.4734188 33.7212879
32 40.6220262 33.4734188
33 37.5906685 40.6220262
34 39.2615968 37.5906685
35 37.9642965 39.2615968
36 39.5581956 37.9642965
37 41.5852061 39.5581956
38 42.9634665 41.5852061
39 40.8465591 42.9634665
40 52.4769361 40.8465591
41 49.4094741 52.4769361
42 42.1798496 49.4094741
43 40.0130546 42.1798496
44 41.2968711 40.0130546
45 24.6977377 41.2968711
46 19.6120596 24.6977377
47 6.1256718 19.6120596
48 17.6923751 6.1256718
49 10.0707924 17.6923751
50 4.0977316 10.0707924
51 2.5483432 4.0977316
52 -1.4427627 2.5483432
53 -15.6454338 -1.4427627
54 -12.6589462 -15.6454338
55 -11.7445815 -12.6589462
56 -29.0822478 -11.7445815
57 -26.4653546 -29.0822478
58 -28.0373920 -26.4653546
59 -29.1723729 -28.0373920
60 -16.7950138 -29.1723729
61 -20.6599045 -16.7950138
62 -29.0920432 -20.6599045
63 -22.1824393 -29.0920432
64 -35.2543626 -22.1824393
65 -28.7543056 -35.2543626
66 -30.4706316 -28.7543056
67 -35.6912194 -30.4706316
68 -43.6777355 -35.6912194
69 -41.0335609 -43.6777355
70 -31.5513778 -41.0335609
71 -20.4703321 -31.5513778
> 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/7lcy91258624273.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/8qbl31258624273.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/95d4q1258624273.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/10cbm51258624273.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/11ct2q1258624273.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/128xko1258624273.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/13kb7a1258624273.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/1414w51258624273.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/158wgt1258624273.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/161p1q1258624273.tab")
+ }
>
> system("convert tmp/11f9e1258624273.ps tmp/11f9e1258624273.png")
> system("convert tmp/2kd3z1258624273.ps tmp/2kd3z1258624273.png")
> system("convert tmp/3wwu11258624273.ps tmp/3wwu11258624273.png")
> system("convert tmp/4idn71258624273.ps tmp/4idn71258624273.png")
> system("convert tmp/5pkv31258624273.ps tmp/5pkv31258624273.png")
> system("convert tmp/6d1b21258624273.ps tmp/6d1b21258624273.png")
> system("convert tmp/7lcy91258624273.ps tmp/7lcy91258624273.png")
> system("convert tmp/8qbl31258624273.ps tmp/8qbl31258624273.png")
> system("convert tmp/95d4q1258624273.ps tmp/95d4q1258624273.png")
> system("convert tmp/10cbm51258624273.ps tmp/10cbm51258624273.png")
>
>
> proc.time()
user system elapsed
2.544 1.565 3.969