R version 2.8.0 (2008-10-20)
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(97.8,0,107.4,0,117.5,0,105.6,0,97.4,0,99.5,0,98.0,0,104.3,0,100.6,0,101.1,0,103.9,0,96.9,0,95.5,0,108.4,0,117.0,0,103.8,0,100.8,0,110.6,0,104.0,0,112.6,0,107.3,0,98.9,0,109.8,0,104.9,0,102.2,0,123.9,0,124.9,0,112.7,0,121.9,0,100.6,0,104.3,0,120.4,0,107.5,0,102.9,0,125.6,0,107.5,0,108.8,0,128.4,0,121.1,0,119.5,0,128.7,0,108.7,0,105.5,0,119.8,0,111.3,0,110.6,0,120.1,0,97.5,0,107.7,0,127.3,0,117.2,0,119.8,0,116.2,0,111.0,0,112.4,0,130.6,0,109.1,0,118.8,0,123.9,0,101.6,0,112.8,0,128.0,0,129.6,0,125.8,0,119.5,0,115.7,0,113.6,0,129.7,0,112.0,0,116.8,0,127.0,0,112.1,1,114.2,1,121.1,1,131.6,1,125.0,1,120.4,1,117.7,1,117.5,1,120.6,1,127.5,1,112.3,1,124.5,1,115.2,1),dim=c(2,84),dimnames=list(c('C','D'),1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('C','D'),1:84))
> 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
C D M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 97.8 0 1 0 0 0 0 0 0 0 0 0 0 1
2 107.4 0 0 1 0 0 0 0 0 0 0 0 0 2
3 117.5 0 0 0 1 0 0 0 0 0 0 0 0 3
4 105.6 0 0 0 0 1 0 0 0 0 0 0 0 4
5 97.4 0 0 0 0 0 1 0 0 0 0 0 0 5
6 99.5 0 0 0 0 0 0 1 0 0 0 0 0 6
7 98.0 0 0 0 0 0 0 0 1 0 0 0 0 7
8 104.3 0 0 0 0 0 0 0 0 1 0 0 0 8
9 100.6 0 0 0 0 0 0 0 0 0 1 0 0 9
10 101.1 0 0 0 0 0 0 0 0 0 0 1 0 10
11 103.9 0 0 0 0 0 0 0 0 0 0 0 1 11
12 96.9 0 0 0 0 0 0 0 0 0 0 0 0 12
13 95.5 0 1 0 0 0 0 0 0 0 0 0 0 13
14 108.4 0 0 1 0 0 0 0 0 0 0 0 0 14
15 117.0 0 0 0 1 0 0 0 0 0 0 0 0 15
16 103.8 0 0 0 0 1 0 0 0 0 0 0 0 16
17 100.8 0 0 0 0 0 1 0 0 0 0 0 0 17
18 110.6 0 0 0 0 0 0 1 0 0 0 0 0 18
19 104.0 0 0 0 0 0 0 0 1 0 0 0 0 19
20 112.6 0 0 0 0 0 0 0 0 1 0 0 0 20
21 107.3 0 0 0 0 0 0 0 0 0 1 0 0 21
22 98.9 0 0 0 0 0 0 0 0 0 0 1 0 22
23 109.8 0 0 0 0 0 0 0 0 0 0 0 1 23
24 104.9 0 0 0 0 0 0 0 0 0 0 0 0 24
25 102.2 0 1 0 0 0 0 0 0 0 0 0 0 25
26 123.9 0 0 1 0 0 0 0 0 0 0 0 0 26
27 124.9 0 0 0 1 0 0 0 0 0 0 0 0 27
28 112.7 0 0 0 0 1 0 0 0 0 0 0 0 28
29 121.9 0 0 0 0 0 1 0 0 0 0 0 0 29
30 100.6 0 0 0 0 0 0 1 0 0 0 0 0 30
31 104.3 0 0 0 0 0 0 0 1 0 0 0 0 31
32 120.4 0 0 0 0 0 0 0 0 1 0 0 0 32
33 107.5 0 0 0 0 0 0 0 0 0 1 0 0 33
34 102.9 0 0 0 0 0 0 0 0 0 0 1 0 34
35 125.6 0 0 0 0 0 0 0 0 0 0 0 1 35
36 107.5 0 0 0 0 0 0 0 0 0 0 0 0 36
37 108.8 0 1 0 0 0 0 0 0 0 0 0 0 37
38 128.4 0 0 1 0 0 0 0 0 0 0 0 0 38
39 121.1 0 0 0 1 0 0 0 0 0 0 0 0 39
40 119.5 0 0 0 0 1 0 0 0 0 0 0 0 40
41 128.7 0 0 0 0 0 1 0 0 0 0 0 0 41
42 108.7 0 0 0 0 0 0 1 0 0 0 0 0 42
43 105.5 0 0 0 0 0 0 0 1 0 0 0 0 43
44 119.8 0 0 0 0 0 0 0 0 1 0 0 0 44
45 111.3 0 0 0 0 0 0 0 0 0 1 0 0 45
46 110.6 0 0 0 0 0 0 0 0 0 0 1 0 46
47 120.1 0 0 0 0 0 0 0 0 0 0 0 1 47
48 97.5 0 0 0 0 0 0 0 0 0 0 0 0 48
49 107.7 0 1 0 0 0 0 0 0 0 0 0 0 49
50 127.3 0 0 1 0 0 0 0 0 0 0 0 0 50
51 117.2 0 0 0 1 0 0 0 0 0 0 0 0 51
52 119.8 0 0 0 0 1 0 0 0 0 0 0 0 52
53 116.2 0 0 0 0 0 1 0 0 0 0 0 0 53
54 111.0 0 0 0 0 0 0 1 0 0 0 0 0 54
55 112.4 0 0 0 0 0 0 0 1 0 0 0 0 55
56 130.6 0 0 0 0 0 0 0 0 1 0 0 0 56
57 109.1 0 0 0 0 0 0 0 0 0 1 0 0 57
58 118.8 0 0 0 0 0 0 0 0 0 0 1 0 58
59 123.9 0 0 0 0 0 0 0 0 0 0 0 1 59
60 101.6 0 0 0 0 0 0 0 0 0 0 0 0 60
61 112.8 0 1 0 0 0 0 0 0 0 0 0 0 61
62 128.0 0 0 1 0 0 0 0 0 0 0 0 0 62
63 129.6 0 0 0 1 0 0 0 0 0 0 0 0 63
64 125.8 0 0 0 0 1 0 0 0 0 0 0 0 64
65 119.5 0 0 0 0 0 1 0 0 0 0 0 0 65
66 115.7 0 0 0 0 0 0 1 0 0 0 0 0 66
67 113.6 0 0 0 0 0 0 0 1 0 0 0 0 67
68 129.7 0 0 0 0 0 0 0 0 1 0 0 0 68
69 112.0 0 0 0 0 0 0 0 0 0 1 0 0 69
70 116.8 0 0 0 0 0 0 0 0 0 0 1 0 70
71 127.0 0 0 0 0 0 0 0 0 0 0 0 1 71
72 112.1 1 0 0 0 0 0 0 0 0 0 0 0 72
73 114.2 1 1 0 0 0 0 0 0 0 0 0 0 73
74 121.1 1 0 1 0 0 0 0 0 0 0 0 0 74
75 131.6 1 0 0 1 0 0 0 0 0 0 0 0 75
76 125.0 1 0 0 0 1 0 0 0 0 0 0 0 76
77 120.4 1 0 0 0 0 1 0 0 0 0 0 0 77
78 117.7 1 0 0 0 0 0 1 0 0 0 0 0 78
79 117.5 1 0 0 0 0 0 0 1 0 0 0 0 79
80 120.6 1 0 0 0 0 0 0 0 1 0 0 0 80
81 127.5 1 0 0 0 0 0 0 0 0 1 0 0 81
82 112.3 1 0 0 0 0 0 0 0 0 0 1 0 82
83 124.5 1 0 0 0 0 0 0 0 0 0 0 1 83
84 115.2 1 0 0 0 0 0 0 0 0 0 0 0 84
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) D M1 M2 M3 M4
92.2386 -3.4541 3.1516 17.9345 19.7031 12.7432
M5 M6 M7 M8 M9 M10
11.4118 5.2519 3.7491 15.2749 6.0292 3.7550
M11 t
13.9522 0.2885
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-9.4555 -2.8958 -0.1223 2.3981 13.2208
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 92.23861 2.16091 42.685 < 2e-16 ***
D -3.45408 1.91207 -1.806 0.0751 .
M1 3.15155 2.63012 1.198 0.2349
M2 17.93447 2.62872 6.823 2.63e-09 ***
M3 19.70311 2.62762 7.498 1.53e-10 ***
M4 12.74318 2.62684 4.851 7.17e-06 ***
M5 11.41181 2.62637 4.345 4.62e-05 ***
M6 5.25188 2.62622 2.000 0.0494 *
M7 3.74909 2.62637 1.427 0.1579
M8 15.27487 2.62684 5.815 1.66e-07 ***
M9 6.02922 2.62762 2.295 0.0248 *
M10 3.75500 2.62872 1.428 0.1576
M11 13.95221 2.63012 5.305 1.26e-06 ***
t 0.28851 0.02864 10.073 2.96e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.897 on 70 degrees of freedom
Multiple R-squared: 0.7884, Adjusted R-squared: 0.7491
F-statistic: 20.07 on 13 and 70 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.07687056 0.1537411 0.9231294
[2,] 0.41219479 0.8243896 0.5878052
[3,] 0.31914693 0.6382939 0.6808531
[4,] 0.29478250 0.5895650 0.7052175
[5,] 0.21936730 0.4387346 0.7806327
[6,] 0.21131843 0.4226369 0.7886816
[7,] 0.17018344 0.3403669 0.8298166
[8,] 0.14675625 0.2935125 0.8532437
[9,] 0.09646896 0.1929379 0.9035310
[10,] 0.22217978 0.4443596 0.7778202
[11,] 0.16920371 0.3384074 0.8307963
[12,] 0.12287206 0.2457441 0.8771279
[13,] 0.47336353 0.9467271 0.5266365
[14,] 0.69720160 0.6055968 0.3027984
[15,] 0.65184035 0.6963193 0.3481596
[16,] 0.60079799 0.7984040 0.3992020
[17,] 0.54883536 0.9023293 0.4511646
[18,] 0.54365804 0.9126839 0.4563420
[19,] 0.66184144 0.6763171 0.3381586
[20,] 0.62442233 0.7511553 0.3755777
[21,] 0.55098217 0.8980357 0.4490178
[22,] 0.56113545 0.8777291 0.4388645
[23,] 0.58891224 0.8221755 0.4110878
[24,] 0.51909741 0.9618052 0.4809026
[25,] 0.83946798 0.3210640 0.1605320
[26,] 0.81203947 0.3759211 0.1879605
[27,] 0.79850480 0.4029904 0.2014952
[28,] 0.74652924 0.5069415 0.2534708
[29,] 0.69501545 0.6099691 0.3049845
[30,] 0.62593567 0.7481287 0.3740643
[31,] 0.56162160 0.8767568 0.4383784
[32,] 0.71738517 0.5652297 0.2826148
[33,] 0.66509600 0.6698080 0.3349040
[34,] 0.64721068 0.7055786 0.3527893
[35,] 0.79487962 0.4102408 0.2051204
[36,] 0.73853013 0.5229397 0.2614699
[37,] 0.68468762 0.6306248 0.3153124
[38,] 0.62503405 0.7499319 0.3749659
[39,] 0.54144840 0.9171032 0.4585516
[40,] 0.59821626 0.8035675 0.4017837
[41,] 0.63094982 0.7381004 0.3690502
[42,] 0.66333853 0.6733229 0.3366615
[43,] 0.61146947 0.7770611 0.3885305
[44,] 0.70273187 0.5945363 0.2972681
[45,] 0.60728545 0.7854291 0.3927146
[46,] 0.59484974 0.8103005 0.4051503
[47,] 0.48196292 0.9639258 0.5180371
[48,] 0.36924621 0.7384924 0.6307538
[49,] 0.26065627 0.5213125 0.7393437
[50,] 0.16272797 0.3254559 0.8372720
[51,] 0.09540660 0.1908132 0.9045934
> postscript(file="/var/www/html/rcomp/tmp/18aum1227775430.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/28nn91227775430.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/3a55y1227775430.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/4mbny1227775430.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/5djl61227775430.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 = 84
Frequency = 1
1 2 3 4 5 6
2.121337827 -3.350090744 4.692766399 -0.535805030 -7.692947887 0.278480684
7 8 9 10 11 12
-0.007233601 -5.521519316 -0.264376458 2.221337827 -5.464376458 1.199325901
13 14 15 16 17 18
-3.640731138 -5.812159710 0.730697433 -5.797873995 -7.755016852 7.916411719
19 20 21 22 23 24
2.530697433 -0.683588281 2.973554576 -3.440731138 -3.026445424 5.737256935
25 26 27 28 29 30
-0.402800104 6.225771325 5.168628468 -0.359942961 9.882914182 -5.545657247
31 32 33 34 35 36
-0.631371532 3.654342753 -0.288514389 -2.902800104 9.311485611 4.875187970
37 38 39 40 41 42
2.735130931 7.263702359 -2.093440498 2.977988074 13.220845216 -0.907726212
43 44 45 46 47 48
-2.893440498 -0.407726212 0.049416645 1.335130931 0.349416645 -8.586880996
49 50 51 52 53 54
-1.826938035 2.701633394 -9.455509463 -0.184080892 -2.741223749 -2.069795178
55 56 57 58 59 60
0.544490537 6.930204822 -5.612652320 6.073061965 0.687347680 -7.948949961
61 62 63 64 65 66
-0.189007000 -0.060435572 -0.517578429 2.353850143 -2.903292715 -0.831864143
67 68 69 70 71 72
-1.717578429 2.568135857 -6.174721286 0.610993000 0.325278714 2.543064558
73 74 75 76 77 78
1.203007519 -6.968421053 1.474436090 1.545864662 -2.011278195 1.160150376
79 80 81 82 83 84
2.174436090 -6.539849624 9.317293233 -3.896992481 -2.182706767 2.180995592
> postscript(file="/var/www/html/rcomp/tmp/6cykh1227775430.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 = 84
Frequency = 1
lag(myerror, k = 1) myerror
0 2.121337827 NA
1 -3.350090744 2.121337827
2 4.692766399 -3.350090744
3 -0.535805030 4.692766399
4 -7.692947887 -0.535805030
5 0.278480684 -7.692947887
6 -0.007233601 0.278480684
7 -5.521519316 -0.007233601
8 -0.264376458 -5.521519316
9 2.221337827 -0.264376458
10 -5.464376458 2.221337827
11 1.199325901 -5.464376458
12 -3.640731138 1.199325901
13 -5.812159710 -3.640731138
14 0.730697433 -5.812159710
15 -5.797873995 0.730697433
16 -7.755016852 -5.797873995
17 7.916411719 -7.755016852
18 2.530697433 7.916411719
19 -0.683588281 2.530697433
20 2.973554576 -0.683588281
21 -3.440731138 2.973554576
22 -3.026445424 -3.440731138
23 5.737256935 -3.026445424
24 -0.402800104 5.737256935
25 6.225771325 -0.402800104
26 5.168628468 6.225771325
27 -0.359942961 5.168628468
28 9.882914182 -0.359942961
29 -5.545657247 9.882914182
30 -0.631371532 -5.545657247
31 3.654342753 -0.631371532
32 -0.288514389 3.654342753
33 -2.902800104 -0.288514389
34 9.311485611 -2.902800104
35 4.875187970 9.311485611
36 2.735130931 4.875187970
37 7.263702359 2.735130931
38 -2.093440498 7.263702359
39 2.977988074 -2.093440498
40 13.220845216 2.977988074
41 -0.907726212 13.220845216
42 -2.893440498 -0.907726212
43 -0.407726212 -2.893440498
44 0.049416645 -0.407726212
45 1.335130931 0.049416645
46 0.349416645 1.335130931
47 -8.586880996 0.349416645
48 -1.826938035 -8.586880996
49 2.701633394 -1.826938035
50 -9.455509463 2.701633394
51 -0.184080892 -9.455509463
52 -2.741223749 -0.184080892
53 -2.069795178 -2.741223749
54 0.544490537 -2.069795178
55 6.930204822 0.544490537
56 -5.612652320 6.930204822
57 6.073061965 -5.612652320
58 0.687347680 6.073061965
59 -7.948949961 0.687347680
60 -0.189007000 -7.948949961
61 -0.060435572 -0.189007000
62 -0.517578429 -0.060435572
63 2.353850143 -0.517578429
64 -2.903292715 2.353850143
65 -0.831864143 -2.903292715
66 -1.717578429 -0.831864143
67 2.568135857 -1.717578429
68 -6.174721286 2.568135857
69 0.610993000 -6.174721286
70 0.325278714 0.610993000
71 2.543064558 0.325278714
72 1.203007519 2.543064558
73 -6.968421053 1.203007519
74 1.474436090 -6.968421053
75 1.545864662 1.474436090
76 -2.011278195 1.545864662
77 1.160150376 -2.011278195
78 2.174436090 1.160150376
79 -6.539849624 2.174436090
80 9.317293233 -6.539849624
81 -3.896992481 9.317293233
82 -2.182706767 -3.896992481
83 2.180995592 -2.182706767
84 NA 2.180995592
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -3.350090744 2.121337827
[2,] 4.692766399 -3.350090744
[3,] -0.535805030 4.692766399
[4,] -7.692947887 -0.535805030
[5,] 0.278480684 -7.692947887
[6,] -0.007233601 0.278480684
[7,] -5.521519316 -0.007233601
[8,] -0.264376458 -5.521519316
[9,] 2.221337827 -0.264376458
[10,] -5.464376458 2.221337827
[11,] 1.199325901 -5.464376458
[12,] -3.640731138 1.199325901
[13,] -5.812159710 -3.640731138
[14,] 0.730697433 -5.812159710
[15,] -5.797873995 0.730697433
[16,] -7.755016852 -5.797873995
[17,] 7.916411719 -7.755016852
[18,] 2.530697433 7.916411719
[19,] -0.683588281 2.530697433
[20,] 2.973554576 -0.683588281
[21,] -3.440731138 2.973554576
[22,] -3.026445424 -3.440731138
[23,] 5.737256935 -3.026445424
[24,] -0.402800104 5.737256935
[25,] 6.225771325 -0.402800104
[26,] 5.168628468 6.225771325
[27,] -0.359942961 5.168628468
[28,] 9.882914182 -0.359942961
[29,] -5.545657247 9.882914182
[30,] -0.631371532 -5.545657247
[31,] 3.654342753 -0.631371532
[32,] -0.288514389 3.654342753
[33,] -2.902800104 -0.288514389
[34,] 9.311485611 -2.902800104
[35,] 4.875187970 9.311485611
[36,] 2.735130931 4.875187970
[37,] 7.263702359 2.735130931
[38,] -2.093440498 7.263702359
[39,] 2.977988074 -2.093440498
[40,] 13.220845216 2.977988074
[41,] -0.907726212 13.220845216
[42,] -2.893440498 -0.907726212
[43,] -0.407726212 -2.893440498
[44,] 0.049416645 -0.407726212
[45,] 1.335130931 0.049416645
[46,] 0.349416645 1.335130931
[47,] -8.586880996 0.349416645
[48,] -1.826938035 -8.586880996
[49,] 2.701633394 -1.826938035
[50,] -9.455509463 2.701633394
[51,] -0.184080892 -9.455509463
[52,] -2.741223749 -0.184080892
[53,] -2.069795178 -2.741223749
[54,] 0.544490537 -2.069795178
[55,] 6.930204822 0.544490537
[56,] -5.612652320 6.930204822
[57,] 6.073061965 -5.612652320
[58,] 0.687347680 6.073061965
[59,] -7.948949961 0.687347680
[60,] -0.189007000 -7.948949961
[61,] -0.060435572 -0.189007000
[62,] -0.517578429 -0.060435572
[63,] 2.353850143 -0.517578429
[64,] -2.903292715 2.353850143
[65,] -0.831864143 -2.903292715
[66,] -1.717578429 -0.831864143
[67,] 2.568135857 -1.717578429
[68,] -6.174721286 2.568135857
[69,] 0.610993000 -6.174721286
[70,] 0.325278714 0.610993000
[71,] 2.543064558 0.325278714
[72,] 1.203007519 2.543064558
[73,] -6.968421053 1.203007519
[74,] 1.474436090 -6.968421053
[75,] 1.545864662 1.474436090
[76,] -2.011278195 1.545864662
[77,] 1.160150376 -2.011278195
[78,] 2.174436090 1.160150376
[79,] -6.539849624 2.174436090
[80,] 9.317293233 -6.539849624
[81,] -3.896992481 9.317293233
[82,] -2.182706767 -3.896992481
[83,] 2.180995592 -2.182706767
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -3.350090744 2.121337827
2 4.692766399 -3.350090744
3 -0.535805030 4.692766399
4 -7.692947887 -0.535805030
5 0.278480684 -7.692947887
6 -0.007233601 0.278480684
7 -5.521519316 -0.007233601
8 -0.264376458 -5.521519316
9 2.221337827 -0.264376458
10 -5.464376458 2.221337827
11 1.199325901 -5.464376458
12 -3.640731138 1.199325901
13 -5.812159710 -3.640731138
14 0.730697433 -5.812159710
15 -5.797873995 0.730697433
16 -7.755016852 -5.797873995
17 7.916411719 -7.755016852
18 2.530697433 7.916411719
19 -0.683588281 2.530697433
20 2.973554576 -0.683588281
21 -3.440731138 2.973554576
22 -3.026445424 -3.440731138
23 5.737256935 -3.026445424
24 -0.402800104 5.737256935
25 6.225771325 -0.402800104
26 5.168628468 6.225771325
27 -0.359942961 5.168628468
28 9.882914182 -0.359942961
29 -5.545657247 9.882914182
30 -0.631371532 -5.545657247
31 3.654342753 -0.631371532
32 -0.288514389 3.654342753
33 -2.902800104 -0.288514389
34 9.311485611 -2.902800104
35 4.875187970 9.311485611
36 2.735130931 4.875187970
37 7.263702359 2.735130931
38 -2.093440498 7.263702359
39 2.977988074 -2.093440498
40 13.220845216 2.977988074
41 -0.907726212 13.220845216
42 -2.893440498 -0.907726212
43 -0.407726212 -2.893440498
44 0.049416645 -0.407726212
45 1.335130931 0.049416645
46 0.349416645 1.335130931
47 -8.586880996 0.349416645
48 -1.826938035 -8.586880996
49 2.701633394 -1.826938035
50 -9.455509463 2.701633394
51 -0.184080892 -9.455509463
52 -2.741223749 -0.184080892
53 -2.069795178 -2.741223749
54 0.544490537 -2.069795178
55 6.930204822 0.544490537
56 -5.612652320 6.930204822
57 6.073061965 -5.612652320
58 0.687347680 6.073061965
59 -7.948949961 0.687347680
60 -0.189007000 -7.948949961
61 -0.060435572 -0.189007000
62 -0.517578429 -0.060435572
63 2.353850143 -0.517578429
64 -2.903292715 2.353850143
65 -0.831864143 -2.903292715
66 -1.717578429 -0.831864143
67 2.568135857 -1.717578429
68 -6.174721286 2.568135857
69 0.610993000 -6.174721286
70 0.325278714 0.610993000
71 2.543064558 0.325278714
72 1.203007519 2.543064558
73 -6.968421053 1.203007519
74 1.474436090 -6.968421053
75 1.545864662 1.474436090
76 -2.011278195 1.545864662
77 1.160150376 -2.011278195
78 2.174436090 1.160150376
79 -6.539849624 2.174436090
80 9.317293233 -6.539849624
81 -3.896992481 9.317293233
82 -2.182706767 -3.896992481
83 2.180995592 -2.182706767
> 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/7g8fa1227775430.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/838i31227775430.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/91cbm1227775430.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/10q74q1227775430.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/11nye61227775430.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/1265le1227775430.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/13f4m81227775431.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/14t12f1227775431.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/15wgh61227775431.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/163nlr1227775431.tab")
+ }
>
> system("convert tmp/18aum1227775430.ps tmp/18aum1227775430.png")
> system("convert tmp/28nn91227775430.ps tmp/28nn91227775430.png")
> system("convert tmp/3a55y1227775430.ps tmp/3a55y1227775430.png")
> system("convert tmp/4mbny1227775430.ps tmp/4mbny1227775430.png")
> system("convert tmp/5djl61227775430.ps tmp/5djl61227775430.png")
> system("convert tmp/6cykh1227775430.ps tmp/6cykh1227775430.png")
> system("convert tmp/7g8fa1227775430.ps tmp/7g8fa1227775430.png")
> system("convert tmp/838i31227775430.ps tmp/838i31227775430.png")
> system("convert tmp/91cbm1227775430.ps tmp/91cbm1227775430.png")
> system("convert tmp/10q74q1227775430.ps tmp/10q74q1227775430.png")
>
>
> proc.time()
user system elapsed
3.012 1.741 4.522