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(128.7,0,136.9,0,156.9,0,109.1,0,122.3,0,123.9,0,90.9,0,77.9,0,120.3,0,118.9,0,125.5,0,98.9,0,102.9,0,105.9,0,117.6,0,113.6,0,115.9,0,118.9,0,77.6,0,81.2,0,123.1,0,136.6,0,112.1,0,95.1,0,96.3,0,105.7,0,115.8,0,105.7,0,105.7,0,111.1,0,82.4,0,60,0,107.3,0,99.3,0,113.5,0,108.9,0,100.2,0,103.9,0,138.7,0,120.2,0,100.2,0,143.2,0,70.9,0,85.2,0,133,0,136.6,0,117.9,0,106.3,0,122.3,0,125.5,0,148.4,0,126.3,0,99.6,0,140.4,0,80.3,0,92.6,0,138.5,0,110.9,0,119.6,0,105,0,109,0,129.4,0,148.6,0,101.4,0,134.8,0,143.7,0,81.6,0,90.3,0,141.5,0,140.7,0,140.2,0,100.2,0,125.7,0,119.6,0,134.7,0,109,0,116.3,0,146.9,0,97.4,0,89.4,0,132.1,1,139.8,1,129,1,112.5,1,121.9,1,121.7,1,123.1,1,131.6,1,119.3,1,132.5,1,98.3,1,85.1,1,131.7,1,129.3,1,90.7,1,78.6,1,68.9,1,79.1,1,83.5,1,74.1,1,59.7,1,93.3,1,61.3,1,56.6,1,98.5,1),dim=c(2,105),dimnames=list(c('autoprod','crisis'),1:105))
> y <- array(NA,dim=c(2,105),dimnames=list(c('autoprod','crisis'),1:105))
> 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
autoprod crisis M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 128.7 0 1 0 0 0 0 0 0 0 0 0 0 1
2 136.9 0 0 1 0 0 0 0 0 0 0 0 0 2
3 156.9 0 0 0 1 0 0 0 0 0 0 0 0 3
4 109.1 0 0 0 0 1 0 0 0 0 0 0 0 4
5 122.3 0 0 0 0 0 1 0 0 0 0 0 0 5
6 123.9 0 0 0 0 0 0 1 0 0 0 0 0 6
7 90.9 0 0 0 0 0 0 0 1 0 0 0 0 7
8 77.9 0 0 0 0 0 0 0 0 1 0 0 0 8
9 120.3 0 0 0 0 0 0 0 0 0 1 0 0 9
10 118.9 0 0 0 0 0 0 0 0 0 0 1 0 10
11 125.5 0 0 0 0 0 0 0 0 0 0 0 1 11
12 98.9 0 0 0 0 0 0 0 0 0 0 0 0 12
13 102.9 0 1 0 0 0 0 0 0 0 0 0 0 13
14 105.9 0 0 1 0 0 0 0 0 0 0 0 0 14
15 117.6 0 0 0 1 0 0 0 0 0 0 0 0 15
16 113.6 0 0 0 0 1 0 0 0 0 0 0 0 16
17 115.9 0 0 0 0 0 1 0 0 0 0 0 0 17
18 118.9 0 0 0 0 0 0 1 0 0 0 0 0 18
19 77.6 0 0 0 0 0 0 0 1 0 0 0 0 19
20 81.2 0 0 0 0 0 0 0 0 1 0 0 0 20
21 123.1 0 0 0 0 0 0 0 0 0 1 0 0 21
22 136.6 0 0 0 0 0 0 0 0 0 0 1 0 22
23 112.1 0 0 0 0 0 0 0 0 0 0 0 1 23
24 95.1 0 0 0 0 0 0 0 0 0 0 0 0 24
25 96.3 0 1 0 0 0 0 0 0 0 0 0 0 25
26 105.7 0 0 1 0 0 0 0 0 0 0 0 0 26
27 115.8 0 0 0 1 0 0 0 0 0 0 0 0 27
28 105.7 0 0 0 0 1 0 0 0 0 0 0 0 28
29 105.7 0 0 0 0 0 1 0 0 0 0 0 0 29
30 111.1 0 0 0 0 0 0 1 0 0 0 0 0 30
31 82.4 0 0 0 0 0 0 0 1 0 0 0 0 31
32 60.0 0 0 0 0 0 0 0 0 1 0 0 0 32
33 107.3 0 0 0 0 0 0 0 0 0 1 0 0 33
34 99.3 0 0 0 0 0 0 0 0 0 0 1 0 34
35 113.5 0 0 0 0 0 0 0 0 0 0 0 1 35
36 108.9 0 0 0 0 0 0 0 0 0 0 0 0 36
37 100.2 0 1 0 0 0 0 0 0 0 0 0 0 37
38 103.9 0 0 1 0 0 0 0 0 0 0 0 0 38
39 138.7 0 0 0 1 0 0 0 0 0 0 0 0 39
40 120.2 0 0 0 0 1 0 0 0 0 0 0 0 40
41 100.2 0 0 0 0 0 1 0 0 0 0 0 0 41
42 143.2 0 0 0 0 0 0 1 0 0 0 0 0 42
43 70.9 0 0 0 0 0 0 0 1 0 0 0 0 43
44 85.2 0 0 0 0 0 0 0 0 1 0 0 0 44
45 133.0 0 0 0 0 0 0 0 0 0 1 0 0 45
46 136.6 0 0 0 0 0 0 0 0 0 0 1 0 46
47 117.9 0 0 0 0 0 0 0 0 0 0 0 1 47
48 106.3 0 0 0 0 0 0 0 0 0 0 0 0 48
49 122.3 0 1 0 0 0 0 0 0 0 0 0 0 49
50 125.5 0 0 1 0 0 0 0 0 0 0 0 0 50
51 148.4 0 0 0 1 0 0 0 0 0 0 0 0 51
52 126.3 0 0 0 0 1 0 0 0 0 0 0 0 52
53 99.6 0 0 0 0 0 1 0 0 0 0 0 0 53
54 140.4 0 0 0 0 0 0 1 0 0 0 0 0 54
55 80.3 0 0 0 0 0 0 0 1 0 0 0 0 55
56 92.6 0 0 0 0 0 0 0 0 1 0 0 0 56
57 138.5 0 0 0 0 0 0 0 0 0 1 0 0 57
58 110.9 0 0 0 0 0 0 0 0 0 0 1 0 58
59 119.6 0 0 0 0 0 0 0 0 0 0 0 1 59
60 105.0 0 0 0 0 0 0 0 0 0 0 0 0 60
61 109.0 0 1 0 0 0 0 0 0 0 0 0 0 61
62 129.4 0 0 1 0 0 0 0 0 0 0 0 0 62
63 148.6 0 0 0 1 0 0 0 0 0 0 0 0 63
64 101.4 0 0 0 0 1 0 0 0 0 0 0 0 64
65 134.8 0 0 0 0 0 1 0 0 0 0 0 0 65
66 143.7 0 0 0 0 0 0 1 0 0 0 0 0 66
67 81.6 0 0 0 0 0 0 0 1 0 0 0 0 67
68 90.3 0 0 0 0 0 0 0 0 1 0 0 0 68
69 141.5 0 0 0 0 0 0 0 0 0 1 0 0 69
70 140.7 0 0 0 0 0 0 0 0 0 0 1 0 70
71 140.2 0 0 0 0 0 0 0 0 0 0 0 1 71
72 100.2 0 0 0 0 0 0 0 0 0 0 0 0 72
73 125.7 0 1 0 0 0 0 0 0 0 0 0 0 73
74 119.6 0 0 1 0 0 0 0 0 0 0 0 0 74
75 134.7 0 0 0 1 0 0 0 0 0 0 0 0 75
76 109.0 0 0 0 0 1 0 0 0 0 0 0 0 76
77 116.3 0 0 0 0 0 1 0 0 0 0 0 0 77
78 146.9 0 0 0 0 0 0 1 0 0 0 0 0 78
79 97.4 0 0 0 0 0 0 0 1 0 0 0 0 79
80 89.4 0 0 0 0 0 0 0 0 1 0 0 0 80
81 132.1 1 0 0 0 0 0 0 0 0 1 0 0 81
82 139.8 1 0 0 0 0 0 0 0 0 0 1 0 82
83 129.0 1 0 0 0 0 0 0 0 0 0 0 1 83
84 112.5 1 0 0 0 0 0 0 0 0 0 0 0 84
85 121.9 1 1 0 0 0 0 0 0 0 0 0 0 85
86 121.7 1 0 1 0 0 0 0 0 0 0 0 0 86
87 123.1 1 0 0 1 0 0 0 0 0 0 0 0 87
88 131.6 1 0 0 0 1 0 0 0 0 0 0 0 88
89 119.3 1 0 0 0 0 1 0 0 0 0 0 0 89
90 132.5 1 0 0 0 0 0 1 0 0 0 0 0 90
91 98.3 1 0 0 0 0 0 0 1 0 0 0 0 91
92 85.1 1 0 0 0 0 0 0 0 1 0 0 0 92
93 131.7 1 0 0 0 0 0 0 0 0 1 0 0 93
94 129.3 1 0 0 0 0 0 0 0 0 0 1 0 94
95 90.7 1 0 0 0 0 0 0 0 0 0 0 1 95
96 78.6 1 0 0 0 0 0 0 0 0 0 0 0 96
97 68.9 1 1 0 0 0 0 0 0 0 0 0 0 97
98 79.1 1 0 1 0 0 0 0 0 0 0 0 0 98
99 83.5 1 0 0 1 0 0 0 0 0 0 0 0 99
100 74.1 1 0 0 0 1 0 0 0 0 0 0 0 100
101 59.7 1 0 0 0 0 1 0 0 0 0 0 0 101
102 93.3 1 0 0 0 0 0 1 0 0 0 0 0 102
103 61.3 1 0 0 0 0 0 0 1 0 0 0 0 103
104 56.6 1 0 0 0 0 0 0 0 1 0 0 0 104
105 98.5 1 0 0 0 0 0 0 0 0 1 0 0 105
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) crisis M1 M2 M3 M4
100.93666 -16.18422 7.64784 13.33308 28.77388 9.11468
M5 M6 M7 M8 M9 M10
7.13325 27.07405 -18.90737 -21.46658 25.56136 25.96563
M11 t
17.94531 0.07031
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-39.2873 -10.3763 0.1271 10.7709 31.5454
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 100.93666 6.45625 15.634 < 2e-16 ***
crisis -16.18422 5.39299 -3.001 0.003473 **
M1 7.64784 7.69049 0.994 0.322639
M2 13.33308 7.68822 1.734 0.086265 .
M3 28.77388 7.68669 3.743 0.000318 ***
M4 9.11468 7.68592 1.186 0.238752
M5 7.13325 7.68589 0.928 0.355812
M6 27.07405 7.68660 3.522 0.000671 ***
M7 -18.90737 7.68807 -2.459 0.015811 *
M8 -21.46658 7.69028 -2.791 0.006395 **
M9 25.56136 7.69182 3.323 0.001283 **
M10 25.96563 7.90940 3.283 0.001459 **
M11 17.94531 7.90831 2.269 0.025621 *
t 0.07031 0.07580 0.928 0.356045
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15.82 on 91 degrees of freedom
Multiple R-squared: 0.5718, Adjusted R-squared: 0.5107
F-statistic: 9.349 on 13 and 91 DF, p-value: 5.83e-12
> 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.511977937 0.976044126 0.48802206
[2,] 0.398160044 0.796320088 0.60183996
[3,] 0.261021748 0.522043497 0.73897825
[4,] 0.235721147 0.471442294 0.76427885
[5,] 0.194157894 0.388315789 0.80584211
[6,] 0.259566465 0.519132930 0.74043354
[7,] 0.179346375 0.358692750 0.82065363
[8,] 0.120929620 0.241859240 0.87907038
[9,] 0.080231940 0.160463880 0.91976806
[10,] 0.049024936 0.098049872 0.95097506
[11,] 0.031601317 0.063202635 0.96839868
[12,] 0.021598354 0.043196709 0.97840165
[13,] 0.012039039 0.024078078 0.98796096
[14,] 0.007727555 0.015455110 0.99227245
[15,] 0.005612517 0.011225034 0.99438748
[16,] 0.003987483 0.007974965 0.99601252
[17,] 0.002772997 0.005545994 0.99722700
[18,] 0.004188827 0.008377653 0.99581117
[19,] 0.003021806 0.006043611 0.99697819
[20,] 0.006415420 0.012830840 0.99358458
[21,] 0.004793852 0.009587704 0.99520615
[22,] 0.003499529 0.006999059 0.99650047
[23,] 0.006049838 0.012099675 0.99395016
[24,] 0.009242788 0.018485576 0.99075721
[25,] 0.006784557 0.013569113 0.99321544
[26,] 0.025313552 0.050627104 0.97468645
[27,] 0.024330031 0.048660062 0.97566997
[28,] 0.029714071 0.059428142 0.97028593
[29,] 0.039354248 0.078708496 0.96064575
[30,] 0.049075491 0.098150983 0.95092451
[31,] 0.042668352 0.085336705 0.95733165
[32,] 0.034962280 0.069924561 0.96503772
[33,] 0.035386933 0.070773867 0.96461307
[34,] 0.031600772 0.063201543 0.96839923
[35,] 0.029253801 0.058507602 0.97074620
[36,] 0.024371756 0.048743512 0.97562824
[37,] 0.029758289 0.059516578 0.97024171
[38,] 0.028886061 0.057772121 0.97111394
[39,] 0.035105424 0.070210848 0.96489458
[40,] 0.037342512 0.074685025 0.96265749
[41,] 0.038394623 0.076789247 0.96160538
[42,] 0.118154256 0.236308513 0.88184574
[43,] 0.134632021 0.269264042 0.86536798
[44,] 0.131585459 0.263170918 0.86841454
[45,] 0.153476217 0.306952433 0.84652378
[46,] 0.135367850 0.270735699 0.86463215
[47,] 0.110328635 0.220657270 0.88967137
[48,] 0.224016513 0.448033027 0.77598349
[49,] 0.227147365 0.454294730 0.77285264
[50,] 0.227706659 0.455413318 0.77229334
[51,] 0.516643423 0.966713153 0.48335658
[52,] 0.683190860 0.633618280 0.31680914
[53,] 0.697917225 0.604165550 0.30208278
[54,] 0.708454403 0.583091194 0.29154560
[55,] 0.687951167 0.624097665 0.31204883
[56,] 0.700445028 0.599109943 0.29955497
[57,] 0.654296919 0.691406162 0.34570308
[58,] 0.583332474 0.833335052 0.41666753
[59,] 0.539499060 0.921001880 0.46050094
[60,] 0.544582627 0.910834745 0.45541737
[61,] 0.463633598 0.927267195 0.53636640
[62,] 0.441356813 0.882713625 0.55864319
[63,] 0.365111475 0.730222950 0.63488853
[64,] 0.283776129 0.567552258 0.71622387
[65,] 0.720490034 0.559019931 0.27950997
[66,] 0.919914756 0.160170488 0.08008524
[67,] 0.871447034 0.257105931 0.12855297
[68,] 0.827600469 0.344799061 0.17239953
[69,] 0.782044517 0.435910966 0.21795548
[70,] 0.666033730 0.667932539 0.33396627
[71,] 0.543593867 0.912812266 0.45640613
[72,] 0.553140147 0.893719705 0.44685985
> postscript(file="/var/www/html/rcomp/tmp/1rp791261332533.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/2t9251261332533.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/3ydbt1261332533.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/42g6v1261332533.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/50oz51261332533.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 = 105
Frequency = 1
1 2 3 4 5 6
20.0451893 22.4896338 26.9785227 -1.2325885 13.8785227 -4.5325885
7 8 9 10 11 12
8.3785227 -2.1325885 -6.8308349 -8.7054159 5.8445841 -2.8804159
13 14 15 16 17 18
-6.5985646 -9.3541201 -13.1652312 2.4236577 6.6347688 -10.3763423
19 20 21 22 23 24
-5.7652312 0.3236577 -4.8745888 8.1508302 -8.3991698 -7.5241698
25 26 27 28 29 30
-14.0423184 -10.3978740 -15.8089851 -6.3200962 -4.4089851 -19.0200962
31 32 33 34 35 36
-1.8089851 -21.7200962 -21.5183426 -29.9929236 -7.8429236 5.4320764
37 38 39 40 41 42
-10.9860723 -13.0416279 6.2472610 7.3361499 -10.7527390 12.2361499
43 44 45 46 47 48
-14.1527390 2.6361499 3.3379035 6.4633225 -4.2866775 1.9883225
49 50 51 52 53 54
10.2701738 7.7146182 15.1035071 12.5923960 -12.1964929 8.5923960
55 56 57 58 59 60
-5.5964929 9.1923960 7.9941496 -20.0804314 -3.4304314 -0.1554314
61 62 63 64 65 66
-3.8735801 10.7708644 14.4597533 -13.1513579 22.1597533 11.0486421
67 68 69 70 71 72
-5.1402467 6.0486421 10.1503957 8.8758147 16.3258147 -5.7991853
73 74 75 76 77 78
11.9826660 0.1271105 -0.2840006 -6.3951117 2.8159994 13.4048883
79 80 81 82 83 84
9.8159994 4.3048883 16.0908597 23.3162787 20.4662787 21.8412787
85 86 87 88 89 90
23.5231300 17.5675745 3.4564634 31.5453523 21.1564634 14.3453523
91 92 93 94 95 96
26.0564634 15.3453523 14.8471058 11.9725248 -18.6774752 -12.9024752
97 98 99 100 101 102
-30.3206238 -25.8761794 -36.9872905 -26.7984016 -39.2872905 -25.6984016
103 104 105
-11.7872905 -13.9984016 -19.1966480
> postscript(file="/var/www/html/rcomp/tmp/6wr7z1261332533.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 = 105
Frequency = 1
lag(myerror, k = 1) myerror
0 20.0451893 NA
1 22.4896338 20.0451893
2 26.9785227 22.4896338
3 -1.2325885 26.9785227
4 13.8785227 -1.2325885
5 -4.5325885 13.8785227
6 8.3785227 -4.5325885
7 -2.1325885 8.3785227
8 -6.8308349 -2.1325885
9 -8.7054159 -6.8308349
10 5.8445841 -8.7054159
11 -2.8804159 5.8445841
12 -6.5985646 -2.8804159
13 -9.3541201 -6.5985646
14 -13.1652312 -9.3541201
15 2.4236577 -13.1652312
16 6.6347688 2.4236577
17 -10.3763423 6.6347688
18 -5.7652312 -10.3763423
19 0.3236577 -5.7652312
20 -4.8745888 0.3236577
21 8.1508302 -4.8745888
22 -8.3991698 8.1508302
23 -7.5241698 -8.3991698
24 -14.0423184 -7.5241698
25 -10.3978740 -14.0423184
26 -15.8089851 -10.3978740
27 -6.3200962 -15.8089851
28 -4.4089851 -6.3200962
29 -19.0200962 -4.4089851
30 -1.8089851 -19.0200962
31 -21.7200962 -1.8089851
32 -21.5183426 -21.7200962
33 -29.9929236 -21.5183426
34 -7.8429236 -29.9929236
35 5.4320764 -7.8429236
36 -10.9860723 5.4320764
37 -13.0416279 -10.9860723
38 6.2472610 -13.0416279
39 7.3361499 6.2472610
40 -10.7527390 7.3361499
41 12.2361499 -10.7527390
42 -14.1527390 12.2361499
43 2.6361499 -14.1527390
44 3.3379035 2.6361499
45 6.4633225 3.3379035
46 -4.2866775 6.4633225
47 1.9883225 -4.2866775
48 10.2701738 1.9883225
49 7.7146182 10.2701738
50 15.1035071 7.7146182
51 12.5923960 15.1035071
52 -12.1964929 12.5923960
53 8.5923960 -12.1964929
54 -5.5964929 8.5923960
55 9.1923960 -5.5964929
56 7.9941496 9.1923960
57 -20.0804314 7.9941496
58 -3.4304314 -20.0804314
59 -0.1554314 -3.4304314
60 -3.8735801 -0.1554314
61 10.7708644 -3.8735801
62 14.4597533 10.7708644
63 -13.1513579 14.4597533
64 22.1597533 -13.1513579
65 11.0486421 22.1597533
66 -5.1402467 11.0486421
67 6.0486421 -5.1402467
68 10.1503957 6.0486421
69 8.8758147 10.1503957
70 16.3258147 8.8758147
71 -5.7991853 16.3258147
72 11.9826660 -5.7991853
73 0.1271105 11.9826660
74 -0.2840006 0.1271105
75 -6.3951117 -0.2840006
76 2.8159994 -6.3951117
77 13.4048883 2.8159994
78 9.8159994 13.4048883
79 4.3048883 9.8159994
80 16.0908597 4.3048883
81 23.3162787 16.0908597
82 20.4662787 23.3162787
83 21.8412787 20.4662787
84 23.5231300 21.8412787
85 17.5675745 23.5231300
86 3.4564634 17.5675745
87 31.5453523 3.4564634
88 21.1564634 31.5453523
89 14.3453523 21.1564634
90 26.0564634 14.3453523
91 15.3453523 26.0564634
92 14.8471058 15.3453523
93 11.9725248 14.8471058
94 -18.6774752 11.9725248
95 -12.9024752 -18.6774752
96 -30.3206238 -12.9024752
97 -25.8761794 -30.3206238
98 -36.9872905 -25.8761794
99 -26.7984016 -36.9872905
100 -39.2872905 -26.7984016
101 -25.6984016 -39.2872905
102 -11.7872905 -25.6984016
103 -13.9984016 -11.7872905
104 -19.1966480 -13.9984016
105 NA -19.1966480
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 22.4896338 20.0451893
[2,] 26.9785227 22.4896338
[3,] -1.2325885 26.9785227
[4,] 13.8785227 -1.2325885
[5,] -4.5325885 13.8785227
[6,] 8.3785227 -4.5325885
[7,] -2.1325885 8.3785227
[8,] -6.8308349 -2.1325885
[9,] -8.7054159 -6.8308349
[10,] 5.8445841 -8.7054159
[11,] -2.8804159 5.8445841
[12,] -6.5985646 -2.8804159
[13,] -9.3541201 -6.5985646
[14,] -13.1652312 -9.3541201
[15,] 2.4236577 -13.1652312
[16,] 6.6347688 2.4236577
[17,] -10.3763423 6.6347688
[18,] -5.7652312 -10.3763423
[19,] 0.3236577 -5.7652312
[20,] -4.8745888 0.3236577
[21,] 8.1508302 -4.8745888
[22,] -8.3991698 8.1508302
[23,] -7.5241698 -8.3991698
[24,] -14.0423184 -7.5241698
[25,] -10.3978740 -14.0423184
[26,] -15.8089851 -10.3978740
[27,] -6.3200962 -15.8089851
[28,] -4.4089851 -6.3200962
[29,] -19.0200962 -4.4089851
[30,] -1.8089851 -19.0200962
[31,] -21.7200962 -1.8089851
[32,] -21.5183426 -21.7200962
[33,] -29.9929236 -21.5183426
[34,] -7.8429236 -29.9929236
[35,] 5.4320764 -7.8429236
[36,] -10.9860723 5.4320764
[37,] -13.0416279 -10.9860723
[38,] 6.2472610 -13.0416279
[39,] 7.3361499 6.2472610
[40,] -10.7527390 7.3361499
[41,] 12.2361499 -10.7527390
[42,] -14.1527390 12.2361499
[43,] 2.6361499 -14.1527390
[44,] 3.3379035 2.6361499
[45,] 6.4633225 3.3379035
[46,] -4.2866775 6.4633225
[47,] 1.9883225 -4.2866775
[48,] 10.2701738 1.9883225
[49,] 7.7146182 10.2701738
[50,] 15.1035071 7.7146182
[51,] 12.5923960 15.1035071
[52,] -12.1964929 12.5923960
[53,] 8.5923960 -12.1964929
[54,] -5.5964929 8.5923960
[55,] 9.1923960 -5.5964929
[56,] 7.9941496 9.1923960
[57,] -20.0804314 7.9941496
[58,] -3.4304314 -20.0804314
[59,] -0.1554314 -3.4304314
[60,] -3.8735801 -0.1554314
[61,] 10.7708644 -3.8735801
[62,] 14.4597533 10.7708644
[63,] -13.1513579 14.4597533
[64,] 22.1597533 -13.1513579
[65,] 11.0486421 22.1597533
[66,] -5.1402467 11.0486421
[67,] 6.0486421 -5.1402467
[68,] 10.1503957 6.0486421
[69,] 8.8758147 10.1503957
[70,] 16.3258147 8.8758147
[71,] -5.7991853 16.3258147
[72,] 11.9826660 -5.7991853
[73,] 0.1271105 11.9826660
[74,] -0.2840006 0.1271105
[75,] -6.3951117 -0.2840006
[76,] 2.8159994 -6.3951117
[77,] 13.4048883 2.8159994
[78,] 9.8159994 13.4048883
[79,] 4.3048883 9.8159994
[80,] 16.0908597 4.3048883
[81,] 23.3162787 16.0908597
[82,] 20.4662787 23.3162787
[83,] 21.8412787 20.4662787
[84,] 23.5231300 21.8412787
[85,] 17.5675745 23.5231300
[86,] 3.4564634 17.5675745
[87,] 31.5453523 3.4564634
[88,] 21.1564634 31.5453523
[89,] 14.3453523 21.1564634
[90,] 26.0564634 14.3453523
[91,] 15.3453523 26.0564634
[92,] 14.8471058 15.3453523
[93,] 11.9725248 14.8471058
[94,] -18.6774752 11.9725248
[95,] -12.9024752 -18.6774752
[96,] -30.3206238 -12.9024752
[97,] -25.8761794 -30.3206238
[98,] -36.9872905 -25.8761794
[99,] -26.7984016 -36.9872905
[100,] -39.2872905 -26.7984016
[101,] -25.6984016 -39.2872905
[102,] -11.7872905 -25.6984016
[103,] -13.9984016 -11.7872905
[104,] -19.1966480 -13.9984016
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 22.4896338 20.0451893
2 26.9785227 22.4896338
3 -1.2325885 26.9785227
4 13.8785227 -1.2325885
5 -4.5325885 13.8785227
6 8.3785227 -4.5325885
7 -2.1325885 8.3785227
8 -6.8308349 -2.1325885
9 -8.7054159 -6.8308349
10 5.8445841 -8.7054159
11 -2.8804159 5.8445841
12 -6.5985646 -2.8804159
13 -9.3541201 -6.5985646
14 -13.1652312 -9.3541201
15 2.4236577 -13.1652312
16 6.6347688 2.4236577
17 -10.3763423 6.6347688
18 -5.7652312 -10.3763423
19 0.3236577 -5.7652312
20 -4.8745888 0.3236577
21 8.1508302 -4.8745888
22 -8.3991698 8.1508302
23 -7.5241698 -8.3991698
24 -14.0423184 -7.5241698
25 -10.3978740 -14.0423184
26 -15.8089851 -10.3978740
27 -6.3200962 -15.8089851
28 -4.4089851 -6.3200962
29 -19.0200962 -4.4089851
30 -1.8089851 -19.0200962
31 -21.7200962 -1.8089851
32 -21.5183426 -21.7200962
33 -29.9929236 -21.5183426
34 -7.8429236 -29.9929236
35 5.4320764 -7.8429236
36 -10.9860723 5.4320764
37 -13.0416279 -10.9860723
38 6.2472610 -13.0416279
39 7.3361499 6.2472610
40 -10.7527390 7.3361499
41 12.2361499 -10.7527390
42 -14.1527390 12.2361499
43 2.6361499 -14.1527390
44 3.3379035 2.6361499
45 6.4633225 3.3379035
46 -4.2866775 6.4633225
47 1.9883225 -4.2866775
48 10.2701738 1.9883225
49 7.7146182 10.2701738
50 15.1035071 7.7146182
51 12.5923960 15.1035071
52 -12.1964929 12.5923960
53 8.5923960 -12.1964929
54 -5.5964929 8.5923960
55 9.1923960 -5.5964929
56 7.9941496 9.1923960
57 -20.0804314 7.9941496
58 -3.4304314 -20.0804314
59 -0.1554314 -3.4304314
60 -3.8735801 -0.1554314
61 10.7708644 -3.8735801
62 14.4597533 10.7708644
63 -13.1513579 14.4597533
64 22.1597533 -13.1513579
65 11.0486421 22.1597533
66 -5.1402467 11.0486421
67 6.0486421 -5.1402467
68 10.1503957 6.0486421
69 8.8758147 10.1503957
70 16.3258147 8.8758147
71 -5.7991853 16.3258147
72 11.9826660 -5.7991853
73 0.1271105 11.9826660
74 -0.2840006 0.1271105
75 -6.3951117 -0.2840006
76 2.8159994 -6.3951117
77 13.4048883 2.8159994
78 9.8159994 13.4048883
79 4.3048883 9.8159994
80 16.0908597 4.3048883
81 23.3162787 16.0908597
82 20.4662787 23.3162787
83 21.8412787 20.4662787
84 23.5231300 21.8412787
85 17.5675745 23.5231300
86 3.4564634 17.5675745
87 31.5453523 3.4564634
88 21.1564634 31.5453523
89 14.3453523 21.1564634
90 26.0564634 14.3453523
91 15.3453523 26.0564634
92 14.8471058 15.3453523
93 11.9725248 14.8471058
94 -18.6774752 11.9725248
95 -12.9024752 -18.6774752
96 -30.3206238 -12.9024752
97 -25.8761794 -30.3206238
98 -36.9872905 -25.8761794
99 -26.7984016 -36.9872905
100 -39.2872905 -26.7984016
101 -25.6984016 -39.2872905
102 -11.7872905 -25.6984016
103 -13.9984016 -11.7872905
104 -19.1966480 -13.9984016
> 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/7xel91261332533.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/8l3ot1261332533.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/97ahe1261332533.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/10e4kj1261332533.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/11f6ek1261332533.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/12ylx21261332533.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/13erv81261332533.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/14qhy71261332533.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/15t8rh1261332533.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/16n7zo1261332533.tab")
+ }
>
> try(system("convert tmp/1rp791261332533.ps tmp/1rp791261332533.png",intern=TRUE))
character(0)
> try(system("convert tmp/2t9251261332533.ps tmp/2t9251261332533.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ydbt1261332533.ps tmp/3ydbt1261332533.png",intern=TRUE))
character(0)
> try(system("convert tmp/42g6v1261332533.ps tmp/42g6v1261332533.png",intern=TRUE))
character(0)
> try(system("convert tmp/50oz51261332533.ps tmp/50oz51261332533.png",intern=TRUE))
character(0)
> try(system("convert tmp/6wr7z1261332533.ps tmp/6wr7z1261332533.png",intern=TRUE))
character(0)
> try(system("convert tmp/7xel91261332533.ps tmp/7xel91261332533.png",intern=TRUE))
character(0)
> try(system("convert tmp/8l3ot1261332533.ps tmp/8l3ot1261332533.png",intern=TRUE))
character(0)
> try(system("convert tmp/97ahe1261332533.ps tmp/97ahe1261332533.png",intern=TRUE))
character(0)
> try(system("convert tmp/10e4kj1261332533.ps tmp/10e4kj1261332533.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.073 1.657 3.742