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(106370
+ ,100.3
+ ,109375
+ ,101.9
+ ,116476
+ ,102.1
+ ,123297
+ ,103.2
+ ,114813
+ ,103.7
+ ,117925
+ ,106.2
+ ,126466
+ ,107.7
+ ,131235
+ ,109.9
+ ,120546
+ ,111.7
+ ,123791
+ ,114.9
+ ,129813
+ ,116
+ ,133463
+ ,118.3
+ ,122987
+ ,120.4
+ ,125418
+ ,126
+ ,130199
+ ,128.1
+ ,133016
+ ,130.1
+ ,121454
+ ,130.8
+ ,122044
+ ,133.6
+ ,128313
+ ,134.2
+ ,131556
+ ,135.5
+ ,120027
+ ,136.2
+ ,123001
+ ,139.1
+ ,130111
+ ,139
+ ,132524
+ ,139.6
+ ,123742
+ ,138.7
+ ,124931
+ ,140.9
+ ,133646
+ ,141.3
+ ,136557
+ ,141.8
+ ,127509
+ ,142
+ ,128945
+ ,144.5
+ ,137191
+ ,144.6
+ ,139716
+ ,145.5
+ ,129083
+ ,146.8
+ ,131604
+ ,149.5
+ ,139413
+ ,149.9
+ ,143125
+ ,150.1
+ ,133948
+ ,150.9
+ ,137116
+ ,152.8
+ ,144864
+ ,153.1
+ ,149277
+ ,154
+ ,138796
+ ,154.9
+ ,143258
+ ,156.9
+ ,150034
+ ,158.4
+ ,154708
+ ,159.7
+ ,144888
+ ,160.2
+ ,148762
+ ,163.2
+ ,156500
+ ,163.7
+ ,161088
+ ,164.4
+ ,152772
+ ,163.7
+ ,158011
+ ,165.5
+ ,163318
+ ,165.6
+ ,169969
+ ,166.8
+ ,162269
+ ,167.5
+ ,165765
+ ,170.6
+ ,170600
+ ,170.9
+ ,174681
+ ,172
+ ,166364
+ ,171.8
+ ,170240
+ ,173.9
+ ,176150
+ ,174
+ ,182056
+ ,173.8
+ ,172218
+ ,173.9
+ ,177856
+ ,176
+ ,182253
+ ,176.6
+ ,188090
+ ,178.2
+ ,176863
+ ,179.2
+ ,183273
+ ,181.3
+ ,187969
+ ,181.8
+ ,194650
+ ,182.9
+ ,183036
+ ,183.8
+ ,189516
+ ,186.3
+ ,193805
+ ,187.4
+ ,200499
+ ,189.2
+ ,188142
+ ,189.7
+ ,193732
+ ,191.9
+ ,197126
+ ,192.6
+ ,205140
+ ,193.7
+ ,191751
+ ,194.2
+ ,196700
+ ,197.6
+ ,199784
+ ,199.3
+ ,207360
+ ,201.4
+ ,196101
+ ,203
+ ,200824
+ ,206.3
+ ,205743
+ ,207.1
+ ,212489
+ ,209.8
+ ,200810
+ ,211.1
+ ,203683
+ ,215.3
+ ,207286
+ ,217.4
+ ,210910
+ ,215.5
+ ,194915
+ ,210.9
+ ,217920
+ ,212.6)
+ ,dim=c(2
+ ,90)
+ ,dimnames=list(c('Y'
+ ,'X')
+ ,1:90))
> y <- array(NA,dim=c(2,90),dimnames=list(c('Y','X'),1:90))
> 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 106370 100.3 1 0 0 0 0 0 0 0 0 0 0 1
2 109375 101.9 0 1 0 0 0 0 0 0 0 0 0 2
3 116476 102.1 0 0 1 0 0 0 0 0 0 0 0 3
4 123297 103.2 0 0 0 1 0 0 0 0 0 0 0 4
5 114813 103.7 0 0 0 0 1 0 0 0 0 0 0 5
6 117925 106.2 0 0 0 0 0 1 0 0 0 0 0 6
7 126466 107.7 0 0 0 0 0 0 1 0 0 0 0 7
8 131235 109.9 0 0 0 0 0 0 0 1 0 0 0 8
9 120546 111.7 0 0 0 0 0 0 0 0 1 0 0 9
10 123791 114.9 0 0 0 0 0 0 0 0 0 1 0 10
11 129813 116.0 0 0 0 0 0 0 0 0 0 0 1 11
12 133463 118.3 0 0 0 0 0 0 0 0 0 0 0 12
13 122987 120.4 1 0 0 0 0 0 0 0 0 0 0 13
14 125418 126.0 0 1 0 0 0 0 0 0 0 0 0 14
15 130199 128.1 0 0 1 0 0 0 0 0 0 0 0 15
16 133016 130.1 0 0 0 1 0 0 0 0 0 0 0 16
17 121454 130.8 0 0 0 0 1 0 0 0 0 0 0 17
18 122044 133.6 0 0 0 0 0 1 0 0 0 0 0 18
19 128313 134.2 0 0 0 0 0 0 1 0 0 0 0 19
20 131556 135.5 0 0 0 0 0 0 0 1 0 0 0 20
21 120027 136.2 0 0 0 0 0 0 0 0 1 0 0 21
22 123001 139.1 0 0 0 0 0 0 0 0 0 1 0 22
23 130111 139.0 0 0 0 0 0 0 0 0 0 0 1 23
24 132524 139.6 0 0 0 0 0 0 0 0 0 0 0 24
25 123742 138.7 1 0 0 0 0 0 0 0 0 0 0 25
26 124931 140.9 0 1 0 0 0 0 0 0 0 0 0 26
27 133646 141.3 0 0 1 0 0 0 0 0 0 0 0 27
28 136557 141.8 0 0 0 1 0 0 0 0 0 0 0 28
29 127509 142.0 0 0 0 0 1 0 0 0 0 0 0 29
30 128945 144.5 0 0 0 0 0 1 0 0 0 0 0 30
31 137191 144.6 0 0 0 0 0 0 1 0 0 0 0 31
32 139716 145.5 0 0 0 0 0 0 0 1 0 0 0 32
33 129083 146.8 0 0 0 0 0 0 0 0 1 0 0 33
34 131604 149.5 0 0 0 0 0 0 0 0 0 1 0 34
35 139413 149.9 0 0 0 0 0 0 0 0 0 0 1 35
36 143125 150.1 0 0 0 0 0 0 0 0 0 0 0 36
37 133948 150.9 1 0 0 0 0 0 0 0 0 0 0 37
38 137116 152.8 0 1 0 0 0 0 0 0 0 0 0 38
39 144864 153.1 0 0 1 0 0 0 0 0 0 0 0 39
40 149277 154.0 0 0 0 1 0 0 0 0 0 0 0 40
41 138796 154.9 0 0 0 0 1 0 0 0 0 0 0 41
42 143258 156.9 0 0 0 0 0 1 0 0 0 0 0 42
43 150034 158.4 0 0 0 0 0 0 1 0 0 0 0 43
44 154708 159.7 0 0 0 0 0 0 0 1 0 0 0 44
45 144888 160.2 0 0 0 0 0 0 0 0 1 0 0 45
46 148762 163.2 0 0 0 0 0 0 0 0 0 1 0 46
47 156500 163.7 0 0 0 0 0 0 0 0 0 0 1 47
48 161088 164.4 0 0 0 0 0 0 0 0 0 0 0 48
49 152772 163.7 1 0 0 0 0 0 0 0 0 0 0 49
50 158011 165.5 0 1 0 0 0 0 0 0 0 0 0 50
51 163318 165.6 0 0 1 0 0 0 0 0 0 0 0 51
52 169969 166.8 0 0 0 1 0 0 0 0 0 0 0 52
53 162269 167.5 0 0 0 0 1 0 0 0 0 0 0 53
54 165765 170.6 0 0 0 0 0 1 0 0 0 0 0 54
55 170600 170.9 0 0 0 0 0 0 1 0 0 0 0 55
56 174681 172.0 0 0 0 0 0 0 0 1 0 0 0 56
57 166364 171.8 0 0 0 0 0 0 0 0 1 0 0 57
58 170240 173.9 0 0 0 0 0 0 0 0 0 1 0 58
59 176150 174.0 0 0 0 0 0 0 0 0 0 0 1 59
60 182056 173.8 0 0 0 0 0 0 0 0 0 0 0 60
61 172218 173.9 1 0 0 0 0 0 0 0 0 0 0 61
62 177856 176.0 0 1 0 0 0 0 0 0 0 0 0 62
63 182253 176.6 0 0 1 0 0 0 0 0 0 0 0 63
64 188090 178.2 0 0 0 1 0 0 0 0 0 0 0 64
65 176863 179.2 0 0 0 0 1 0 0 0 0 0 0 65
66 183273 181.3 0 0 0 0 0 1 0 0 0 0 0 66
67 187969 181.8 0 0 0 0 0 0 1 0 0 0 0 67
68 194650 182.9 0 0 0 0 0 0 0 1 0 0 0 68
69 183036 183.8 0 0 0 0 0 0 0 0 1 0 0 69
70 189516 186.3 0 0 0 0 0 0 0 0 0 1 0 70
71 193805 187.4 0 0 0 0 0 0 0 0 0 0 1 71
72 200499 189.2 0 0 0 0 0 0 0 0 0 0 0 72
73 188142 189.7 1 0 0 0 0 0 0 0 0 0 0 73
74 193732 191.9 0 1 0 0 0 0 0 0 0 0 0 74
75 197126 192.6 0 0 1 0 0 0 0 0 0 0 0 75
76 205140 193.7 0 0 0 1 0 0 0 0 0 0 0 76
77 191751 194.2 0 0 0 0 1 0 0 0 0 0 0 77
78 196700 197.6 0 0 0 0 0 1 0 0 0 0 0 78
79 199784 199.3 0 0 0 0 0 0 1 0 0 0 0 79
80 207360 201.4 0 0 0 0 0 0 0 1 0 0 0 80
81 196101 203.0 0 0 0 0 0 0 0 0 1 0 0 81
82 200824 206.3 0 0 0 0 0 0 0 0 0 1 0 82
83 205743 207.1 0 0 0 0 0 0 0 0 0 0 1 83
84 212489 209.8 0 0 0 0 0 0 0 0 0 0 0 84
85 200810 211.1 1 0 0 0 0 0 0 0 0 0 0 85
86 203683 215.3 0 1 0 0 0 0 0 0 0 0 0 86
87 207286 217.4 0 0 1 0 0 0 0 0 0 0 0 87
88 210910 215.5 0 0 0 1 0 0 0 0 0 0 0 88
89 194915 210.9 0 0 0 0 1 0 0 0 0 0 0 89
90 217920 212.6 0 0 0 0 0 1 0 0 0 0 0 90
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
208567.1 -921.0 -11949.2 -8082.8 -3965.7 -343.3
M5 M6 M7 M8 M9 M10
-13602.5 -7617.9 -4314.9 -468.3 -12413.5 -8127.3
M11 t
-3619.4 2261.9
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-8628.8 -3970.5 -108.1 3421.7 9649.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 208567.1 16132.5 12.928 < 2e-16 ***
X -921.0 150.6 -6.115 3.86e-08 ***
M1 -11949.2 2732.8 -4.373 3.85e-05 ***
M2 -8082.8 2722.4 -2.969 0.00400 **
M3 -3965.7 2722.5 -1.457 0.14934
M4 -343.3 2724.0 -0.126 0.90005
M5 -13602.5 2737.8 -4.968 4.05e-06 ***
M6 -7617.9 2723.0 -2.798 0.00652 **
M7 -4314.9 2815.7 -1.532 0.12956
M8 -468.3 2813.4 -0.166 0.86823
M9 -12413.5 2814.5 -4.410 3.35e-05 ***
M10 -8127.3 2812.5 -2.890 0.00502 **
M11 -3619.4 2810.4 -1.288 0.20170
t 2261.9 181.7 12.448 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5258 on 76 degrees of freedom
Multiple R-squared: 0.9757, Adjusted R-squared: 0.9716
F-statistic: 234.9 on 13 and 76 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.2112062 0.422412432 0.788793784
[2,] 0.2610144 0.522028773 0.738985614
[3,] 0.4911433 0.982286657 0.508856672
[4,] 0.7707101 0.458579815 0.229289907
[5,] 0.9374073 0.125185333 0.062592666
[6,] 0.9808135 0.038373011 0.019186505
[7,] 0.9941894 0.011621206 0.005810603
[8,] 0.9975865 0.004826974 0.002413487
[9,] 0.9973710 0.005257969 0.002628985
[10,] 0.9957843 0.008431390 0.004215695
[11,] 0.9964478 0.007104326 0.003552163
[12,] 0.9952537 0.009492670 0.004746335
[13,] 0.9974648 0.005070444 0.002535222
[14,] 0.9954695 0.009060994 0.004530497
[15,] 0.9948713 0.010257434 0.005128717
[16,] 0.9914942 0.017011538 0.008505769
[17,] 0.9862607 0.027478644 0.013739322
[18,] 0.9790646 0.041870831 0.020935416
[19,] 0.9679066 0.064186890 0.032093445
[20,] 0.9582438 0.083512414 0.041756207
[21,] 0.9517114 0.096577150 0.048288575
[22,] 0.9548672 0.090265508 0.045132754
[23,] 0.9486385 0.102723014 0.051361507
[24,] 0.9420465 0.115907091 0.057953546
[25,] 0.9258902 0.148219533 0.074109767
[26,] 0.9471093 0.105781372 0.052890686
[27,] 0.9363785 0.127243036 0.063621518
[28,] 0.9357730 0.128454041 0.064227020
[29,] 0.9370397 0.125920687 0.062960344
[30,] 0.9491505 0.101699050 0.050849525
[31,] 0.9494498 0.101100481 0.050550241
[32,] 0.9660088 0.067982392 0.033991196
[33,] 0.9757380 0.048524057 0.024262029
[34,] 0.9860065 0.027986918 0.013993459
[35,] 0.9858864 0.028227226 0.014113613
[36,] 0.9867210 0.026558067 0.013279034
[37,] 0.9969659 0.006068244 0.003034122
[38,] 0.9974328 0.005134307 0.002567153
[39,] 0.9968151 0.006369783 0.003184891
[40,] 0.9956052 0.008789586 0.004394793
[41,] 0.9939464 0.012107176 0.006053588
[42,] 0.9919303 0.016139446 0.008069723
[43,] 0.9875032 0.024993677 0.012496838
[44,] 0.9873856 0.025228710 0.012614355
[45,] 0.9832794 0.033441149 0.016720575
[46,] 0.9776583 0.044683334 0.022341667
[47,] 0.9660461 0.067907778 0.033953889
[48,] 0.9477959 0.104408114 0.052204057
[49,] 0.9485176 0.102964746 0.051482373
[50,] 0.9342200 0.131560059 0.065780029
[51,] 0.8937169 0.212566213 0.106283107
[52,] 0.8398776 0.320244809 0.160122405
[53,] 0.7651024 0.469795295 0.234897648
[54,] 0.6718396 0.656320793 0.328160396
[55,] 0.5492308 0.901538479 0.450769239
[56,] 0.4305187 0.861037304 0.569481348
[57,] 0.3080035 0.616007099 0.691996451
> postscript(file="/var/www/html/rcomp/tmp/1tbf01258927376.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/2g68f1258927376.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/3a3691258927376.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/4pks61258927376.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/5p0q01258927376.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 = 90
Frequency = 1
1 2 3 4 5 6
-133.08326 -1782.81328 -876.67858 1073.11024 4046.87507 1214.86252
7 8 9 10 11 12
5572.41967 6259.19462 6911.19853 6555.30029 6820.55991 6707.56512
13 14 15 16 17 18
7852.94166 9313.22989 9649.27327 8423.96619 8503.93194 3426.22075
19 20 21 22 23 24
4682.87379 3014.74463 1813.64353 910.44392 1158.49806 -1757.20448
25 26 27 28 29 30
-1680.84164 -4593.96892 -1889.63330 -4402.44723 -2268.98376 -6776.99631
31 32 33 34 35 36
-4003.84555 -6758.37654 -6510.87491 -8051.27543 -6643.71900 -8628.82337
37 38 39 40 41 42
-7381.75277 -8592.18142 -6946.94626 -7589.35836 -6244.19170 -8186.70653
43 44 45 46 47 48
-5594.14938 -5831.27854 -5507.58056 -5418.67972 -3990.02283 -4638.62492
49 50 51 52 53 54
-3912.06116 -3143.59027 -4123.55602 -2251.66675 1690.29900 -205.11082
55 56 57 58 59 60
-658.75915 -1673.08922 -491.09443 -1229.09769 -1996.84264 -2156.34883
61 62 63 64 65 66
-2214.98142 -771.20916 -2200.67263 -774.38154 -83.11442 14.47120
67 68 69 70 71 72
-393.97621 1191.69372 89.79352 2324.19209 856.45171 3326.95464
73 74 75 76 77 78
1117.72387 2605.59659 265.23358 3408.02239 1476.78723 1310.67878
79 80 81 82 83 84
395.43684 3797.11133 3694.91433 4909.11655 3795.07480 7146.48184
85 86 87 88 89 90
6352.05473 6964.93657 6122.97995 2112.75507 -7121.60337 9202.58042
> postscript(file="/var/www/html/rcomp/tmp/6k3vp1258927376.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 = 90
Frequency = 1
lag(myerror, k = 1) myerror
0 -133.08326 NA
1 -1782.81328 -133.08326
2 -876.67858 -1782.81328
3 1073.11024 -876.67858
4 4046.87507 1073.11024
5 1214.86252 4046.87507
6 5572.41967 1214.86252
7 6259.19462 5572.41967
8 6911.19853 6259.19462
9 6555.30029 6911.19853
10 6820.55991 6555.30029
11 6707.56512 6820.55991
12 7852.94166 6707.56512
13 9313.22989 7852.94166
14 9649.27327 9313.22989
15 8423.96619 9649.27327
16 8503.93194 8423.96619
17 3426.22075 8503.93194
18 4682.87379 3426.22075
19 3014.74463 4682.87379
20 1813.64353 3014.74463
21 910.44392 1813.64353
22 1158.49806 910.44392
23 -1757.20448 1158.49806
24 -1680.84164 -1757.20448
25 -4593.96892 -1680.84164
26 -1889.63330 -4593.96892
27 -4402.44723 -1889.63330
28 -2268.98376 -4402.44723
29 -6776.99631 -2268.98376
30 -4003.84555 -6776.99631
31 -6758.37654 -4003.84555
32 -6510.87491 -6758.37654
33 -8051.27543 -6510.87491
34 -6643.71900 -8051.27543
35 -8628.82337 -6643.71900
36 -7381.75277 -8628.82337
37 -8592.18142 -7381.75277
38 -6946.94626 -8592.18142
39 -7589.35836 -6946.94626
40 -6244.19170 -7589.35836
41 -8186.70653 -6244.19170
42 -5594.14938 -8186.70653
43 -5831.27854 -5594.14938
44 -5507.58056 -5831.27854
45 -5418.67972 -5507.58056
46 -3990.02283 -5418.67972
47 -4638.62492 -3990.02283
48 -3912.06116 -4638.62492
49 -3143.59027 -3912.06116
50 -4123.55602 -3143.59027
51 -2251.66675 -4123.55602
52 1690.29900 -2251.66675
53 -205.11082 1690.29900
54 -658.75915 -205.11082
55 -1673.08922 -658.75915
56 -491.09443 -1673.08922
57 -1229.09769 -491.09443
58 -1996.84264 -1229.09769
59 -2156.34883 -1996.84264
60 -2214.98142 -2156.34883
61 -771.20916 -2214.98142
62 -2200.67263 -771.20916
63 -774.38154 -2200.67263
64 -83.11442 -774.38154
65 14.47120 -83.11442
66 -393.97621 14.47120
67 1191.69372 -393.97621
68 89.79352 1191.69372
69 2324.19209 89.79352
70 856.45171 2324.19209
71 3326.95464 856.45171
72 1117.72387 3326.95464
73 2605.59659 1117.72387
74 265.23358 2605.59659
75 3408.02239 265.23358
76 1476.78723 3408.02239
77 1310.67878 1476.78723
78 395.43684 1310.67878
79 3797.11133 395.43684
80 3694.91433 3797.11133
81 4909.11655 3694.91433
82 3795.07480 4909.11655
83 7146.48184 3795.07480
84 6352.05473 7146.48184
85 6964.93657 6352.05473
86 6122.97995 6964.93657
87 2112.75507 6122.97995
88 -7121.60337 2112.75507
89 9202.58042 -7121.60337
90 NA 9202.58042
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1782.81328 -133.08326
[2,] -876.67858 -1782.81328
[3,] 1073.11024 -876.67858
[4,] 4046.87507 1073.11024
[5,] 1214.86252 4046.87507
[6,] 5572.41967 1214.86252
[7,] 6259.19462 5572.41967
[8,] 6911.19853 6259.19462
[9,] 6555.30029 6911.19853
[10,] 6820.55991 6555.30029
[11,] 6707.56512 6820.55991
[12,] 7852.94166 6707.56512
[13,] 9313.22989 7852.94166
[14,] 9649.27327 9313.22989
[15,] 8423.96619 9649.27327
[16,] 8503.93194 8423.96619
[17,] 3426.22075 8503.93194
[18,] 4682.87379 3426.22075
[19,] 3014.74463 4682.87379
[20,] 1813.64353 3014.74463
[21,] 910.44392 1813.64353
[22,] 1158.49806 910.44392
[23,] -1757.20448 1158.49806
[24,] -1680.84164 -1757.20448
[25,] -4593.96892 -1680.84164
[26,] -1889.63330 -4593.96892
[27,] -4402.44723 -1889.63330
[28,] -2268.98376 -4402.44723
[29,] -6776.99631 -2268.98376
[30,] -4003.84555 -6776.99631
[31,] -6758.37654 -4003.84555
[32,] -6510.87491 -6758.37654
[33,] -8051.27543 -6510.87491
[34,] -6643.71900 -8051.27543
[35,] -8628.82337 -6643.71900
[36,] -7381.75277 -8628.82337
[37,] -8592.18142 -7381.75277
[38,] -6946.94626 -8592.18142
[39,] -7589.35836 -6946.94626
[40,] -6244.19170 -7589.35836
[41,] -8186.70653 -6244.19170
[42,] -5594.14938 -8186.70653
[43,] -5831.27854 -5594.14938
[44,] -5507.58056 -5831.27854
[45,] -5418.67972 -5507.58056
[46,] -3990.02283 -5418.67972
[47,] -4638.62492 -3990.02283
[48,] -3912.06116 -4638.62492
[49,] -3143.59027 -3912.06116
[50,] -4123.55602 -3143.59027
[51,] -2251.66675 -4123.55602
[52,] 1690.29900 -2251.66675
[53,] -205.11082 1690.29900
[54,] -658.75915 -205.11082
[55,] -1673.08922 -658.75915
[56,] -491.09443 -1673.08922
[57,] -1229.09769 -491.09443
[58,] -1996.84264 -1229.09769
[59,] -2156.34883 -1996.84264
[60,] -2214.98142 -2156.34883
[61,] -771.20916 -2214.98142
[62,] -2200.67263 -771.20916
[63,] -774.38154 -2200.67263
[64,] -83.11442 -774.38154
[65,] 14.47120 -83.11442
[66,] -393.97621 14.47120
[67,] 1191.69372 -393.97621
[68,] 89.79352 1191.69372
[69,] 2324.19209 89.79352
[70,] 856.45171 2324.19209
[71,] 3326.95464 856.45171
[72,] 1117.72387 3326.95464
[73,] 2605.59659 1117.72387
[74,] 265.23358 2605.59659
[75,] 3408.02239 265.23358
[76,] 1476.78723 3408.02239
[77,] 1310.67878 1476.78723
[78,] 395.43684 1310.67878
[79,] 3797.11133 395.43684
[80,] 3694.91433 3797.11133
[81,] 4909.11655 3694.91433
[82,] 3795.07480 4909.11655
[83,] 7146.48184 3795.07480
[84,] 6352.05473 7146.48184
[85,] 6964.93657 6352.05473
[86,] 6122.97995 6964.93657
[87,] 2112.75507 6122.97995
[88,] -7121.60337 2112.75507
[89,] 9202.58042 -7121.60337
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1782.81328 -133.08326
2 -876.67858 -1782.81328
3 1073.11024 -876.67858
4 4046.87507 1073.11024
5 1214.86252 4046.87507
6 5572.41967 1214.86252
7 6259.19462 5572.41967
8 6911.19853 6259.19462
9 6555.30029 6911.19853
10 6820.55991 6555.30029
11 6707.56512 6820.55991
12 7852.94166 6707.56512
13 9313.22989 7852.94166
14 9649.27327 9313.22989
15 8423.96619 9649.27327
16 8503.93194 8423.96619
17 3426.22075 8503.93194
18 4682.87379 3426.22075
19 3014.74463 4682.87379
20 1813.64353 3014.74463
21 910.44392 1813.64353
22 1158.49806 910.44392
23 -1757.20448 1158.49806
24 -1680.84164 -1757.20448
25 -4593.96892 -1680.84164
26 -1889.63330 -4593.96892
27 -4402.44723 -1889.63330
28 -2268.98376 -4402.44723
29 -6776.99631 -2268.98376
30 -4003.84555 -6776.99631
31 -6758.37654 -4003.84555
32 -6510.87491 -6758.37654
33 -8051.27543 -6510.87491
34 -6643.71900 -8051.27543
35 -8628.82337 -6643.71900
36 -7381.75277 -8628.82337
37 -8592.18142 -7381.75277
38 -6946.94626 -8592.18142
39 -7589.35836 -6946.94626
40 -6244.19170 -7589.35836
41 -8186.70653 -6244.19170
42 -5594.14938 -8186.70653
43 -5831.27854 -5594.14938
44 -5507.58056 -5831.27854
45 -5418.67972 -5507.58056
46 -3990.02283 -5418.67972
47 -4638.62492 -3990.02283
48 -3912.06116 -4638.62492
49 -3143.59027 -3912.06116
50 -4123.55602 -3143.59027
51 -2251.66675 -4123.55602
52 1690.29900 -2251.66675
53 -205.11082 1690.29900
54 -658.75915 -205.11082
55 -1673.08922 -658.75915
56 -491.09443 -1673.08922
57 -1229.09769 -491.09443
58 -1996.84264 -1229.09769
59 -2156.34883 -1996.84264
60 -2214.98142 -2156.34883
61 -771.20916 -2214.98142
62 -2200.67263 -771.20916
63 -774.38154 -2200.67263
64 -83.11442 -774.38154
65 14.47120 -83.11442
66 -393.97621 14.47120
67 1191.69372 -393.97621
68 89.79352 1191.69372
69 2324.19209 89.79352
70 856.45171 2324.19209
71 3326.95464 856.45171
72 1117.72387 3326.95464
73 2605.59659 1117.72387
74 265.23358 2605.59659
75 3408.02239 265.23358
76 1476.78723 3408.02239
77 1310.67878 1476.78723
78 395.43684 1310.67878
79 3797.11133 395.43684
80 3694.91433 3797.11133
81 4909.11655 3694.91433
82 3795.07480 4909.11655
83 7146.48184 3795.07480
84 6352.05473 7146.48184
85 6964.93657 6352.05473
86 6122.97995 6964.93657
87 2112.75507 6122.97995
88 -7121.60337 2112.75507
89 9202.58042 -7121.60337
> 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/748a51258927376.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/8b3cd1258927376.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/9ceqn1258927376.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/107y821258927376.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/11q86j1258927376.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/126q6s1258927376.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/13t3f81258927376.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/1468iv1258927376.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/15yqwh1258927376.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/16ft4b1258927376.tab")
+ }
>
> system("convert tmp/1tbf01258927376.ps tmp/1tbf01258927376.png")
> system("convert tmp/2g68f1258927376.ps tmp/2g68f1258927376.png")
> system("convert tmp/3a3691258927376.ps tmp/3a3691258927376.png")
> system("convert tmp/4pks61258927376.ps tmp/4pks61258927376.png")
> system("convert tmp/5p0q01258927376.ps tmp/5p0q01258927376.png")
> system("convert tmp/6k3vp1258927376.ps tmp/6k3vp1258927376.png")
> system("convert tmp/748a51258927376.ps tmp/748a51258927376.png")
> system("convert tmp/8b3cd1258927376.ps tmp/8b3cd1258927376.png")
> system("convert tmp/9ceqn1258927376.ps tmp/9ceqn1258927376.png")
> system("convert tmp/107y821258927376.ps tmp/107y821258927376.png")
>
>
> proc.time()
user system elapsed
2.793 1.597 10.051