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(111.6,0,104.6,0,91.6,0,98.3,0,97.7,0,106.3,0,102.3,0,106.6,0,108.1,0,93.8,0,88.2,0,108.9,0,114.2,0,102.5,0,94.2,0,97.4,0,98.5,0,106.5,0,102.9,0,97.1,0,103.7,0,93.4,0,85.8,0,108.6,0,110.2,0,101.2,0,101.2,0,96.9,0,99.4,0,118.7,0,108.0,0,101.2,0,119.9,0,94.8,0,95.3,0,118.0,0,115.9,0,111.4,0,108.2,0,108.8,0,109.5,0,124.8,0,115.3,0,109.5,0,124.2,0,92.9,0,98.4,0,120.9,0,111.7,0,116.1,0,109.4,0,111.7,0,114.3,0,133.7,0,114.3,0,126.5,0,131.0,0,104.0,0,108.9,0,128.5,0,132.4,0,128.0,0,116.4,0,120.9,0,118.6,0,133.1,0,121.1,0,127.6,0,135.4,0,114.9,0,114.3,0,128.9,0,138.9,0,129.4,0,115.0,0,128.0,1,127.0,1,128.8,1,137.9,1,128.4,1,135.9,1,122.2,1,113.1,1,136.2,1,138.0,1,115.2,1,111.0,1,99.2,1,102.4,1,112.7,1,105.5,1,98.3,1,116.4,1,97.4,1,93.3,1,117.4,1),dim=c(2,96),dimnames=list(c('Yt','Xt_dummy'),1:96))
> y <- array(NA,dim=c(2,96),dimnames=list(c('Yt','Xt_dummy'),1:96))
> 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
Yt Xt_dummy M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 111.6 0 1 0 0 0 0 0 0 0 0 0 0 1
2 104.6 0 0 1 0 0 0 0 0 0 0 0 0 2
3 91.6 0 0 0 1 0 0 0 0 0 0 0 0 3
4 98.3 0 0 0 0 1 0 0 0 0 0 0 0 4
5 97.7 0 0 0 0 0 1 0 0 0 0 0 0 5
6 106.3 0 0 0 0 0 0 1 0 0 0 0 0 6
7 102.3 0 0 0 0 0 0 0 1 0 0 0 0 7
8 106.6 0 0 0 0 0 0 0 0 1 0 0 0 8
9 108.1 0 0 0 0 0 0 0 0 0 1 0 0 9
10 93.8 0 0 0 0 0 0 0 0 0 0 1 0 10
11 88.2 0 0 0 0 0 0 0 0 0 0 0 1 11
12 108.9 0 0 0 0 0 0 0 0 0 0 0 0 12
13 114.2 0 1 0 0 0 0 0 0 0 0 0 0 13
14 102.5 0 0 1 0 0 0 0 0 0 0 0 0 14
15 94.2 0 0 0 1 0 0 0 0 0 0 0 0 15
16 97.4 0 0 0 0 1 0 0 0 0 0 0 0 16
17 98.5 0 0 0 0 0 1 0 0 0 0 0 0 17
18 106.5 0 0 0 0 0 0 1 0 0 0 0 0 18
19 102.9 0 0 0 0 0 0 0 1 0 0 0 0 19
20 97.1 0 0 0 0 0 0 0 0 1 0 0 0 20
21 103.7 0 0 0 0 0 0 0 0 0 1 0 0 21
22 93.4 0 0 0 0 0 0 0 0 0 0 1 0 22
23 85.8 0 0 0 0 0 0 0 0 0 0 0 1 23
24 108.6 0 0 0 0 0 0 0 0 0 0 0 0 24
25 110.2 0 1 0 0 0 0 0 0 0 0 0 0 25
26 101.2 0 0 1 0 0 0 0 0 0 0 0 0 26
27 101.2 0 0 0 1 0 0 0 0 0 0 0 0 27
28 96.9 0 0 0 0 1 0 0 0 0 0 0 0 28
29 99.4 0 0 0 0 0 1 0 0 0 0 0 0 29
30 118.7 0 0 0 0 0 0 1 0 0 0 0 0 30
31 108.0 0 0 0 0 0 0 0 1 0 0 0 0 31
32 101.2 0 0 0 0 0 0 0 0 1 0 0 0 32
33 119.9 0 0 0 0 0 0 0 0 0 1 0 0 33
34 94.8 0 0 0 0 0 0 0 0 0 0 1 0 34
35 95.3 0 0 0 0 0 0 0 0 0 0 0 1 35
36 118.0 0 0 0 0 0 0 0 0 0 0 0 0 36
37 115.9 0 1 0 0 0 0 0 0 0 0 0 0 37
38 111.4 0 0 1 0 0 0 0 0 0 0 0 0 38
39 108.2 0 0 0 1 0 0 0 0 0 0 0 0 39
40 108.8 0 0 0 0 1 0 0 0 0 0 0 0 40
41 109.5 0 0 0 0 0 1 0 0 0 0 0 0 41
42 124.8 0 0 0 0 0 0 1 0 0 0 0 0 42
43 115.3 0 0 0 0 0 0 0 1 0 0 0 0 43
44 109.5 0 0 0 0 0 0 0 0 1 0 0 0 44
45 124.2 0 0 0 0 0 0 0 0 0 1 0 0 45
46 92.9 0 0 0 0 0 0 0 0 0 0 1 0 46
47 98.4 0 0 0 0 0 0 0 0 0 0 0 1 47
48 120.9 0 0 0 0 0 0 0 0 0 0 0 0 48
49 111.7 0 1 0 0 0 0 0 0 0 0 0 0 49
50 116.1 0 0 1 0 0 0 0 0 0 0 0 0 50
51 109.4 0 0 0 1 0 0 0 0 0 0 0 0 51
52 111.7 0 0 0 0 1 0 0 0 0 0 0 0 52
53 114.3 0 0 0 0 0 1 0 0 0 0 0 0 53
54 133.7 0 0 0 0 0 0 1 0 0 0 0 0 54
55 114.3 0 0 0 0 0 0 0 1 0 0 0 0 55
56 126.5 0 0 0 0 0 0 0 0 1 0 0 0 56
57 131.0 0 0 0 0 0 0 0 0 0 1 0 0 57
58 104.0 0 0 0 0 0 0 0 0 0 0 1 0 58
59 108.9 0 0 0 0 0 0 0 0 0 0 0 1 59
60 128.5 0 0 0 0 0 0 0 0 0 0 0 0 60
61 132.4 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 116.4 0 0 0 1 0 0 0 0 0 0 0 0 63
64 120.9 0 0 0 0 1 0 0 0 0 0 0 0 64
65 118.6 0 0 0 0 0 1 0 0 0 0 0 0 65
66 133.1 0 0 0 0 0 0 1 0 0 0 0 0 66
67 121.1 0 0 0 0 0 0 0 1 0 0 0 0 67
68 127.6 0 0 0 0 0 0 0 0 1 0 0 0 68
69 135.4 0 0 0 0 0 0 0 0 0 1 0 0 69
70 114.9 0 0 0 0 0 0 0 0 0 0 1 0 70
71 114.3 0 0 0 0 0 0 0 0 0 0 0 1 71
72 128.9 0 0 0 0 0 0 0 0 0 0 0 0 72
73 138.9 0 1 0 0 0 0 0 0 0 0 0 0 73
74 129.4 0 0 1 0 0 0 0 0 0 0 0 0 74
75 115.0 0 0 0 1 0 0 0 0 0 0 0 0 75
76 128.0 1 0 0 0 1 0 0 0 0 0 0 0 76
77 127.0 1 0 0 0 0 1 0 0 0 0 0 0 77
78 128.8 1 0 0 0 0 0 1 0 0 0 0 0 78
79 137.9 1 0 0 0 0 0 0 1 0 0 0 0 79
80 128.4 1 0 0 0 0 0 0 0 1 0 0 0 80
81 135.9 1 0 0 0 0 0 0 0 0 1 0 0 81
82 122.2 1 0 0 0 0 0 0 0 0 0 1 0 82
83 113.1 1 0 0 0 0 0 0 0 0 0 0 1 83
84 136.2 1 0 0 0 0 0 0 0 0 0 0 0 84
85 138.0 1 1 0 0 0 0 0 0 0 0 0 0 85
86 115.2 1 0 1 0 0 0 0 0 0 0 0 0 86
87 111.0 1 0 0 1 0 0 0 0 0 0 0 0 87
88 99.2 1 0 0 0 1 0 0 0 0 0 0 0 88
89 102.4 1 0 0 0 0 1 0 0 0 0 0 0 89
90 112.7 1 0 0 0 0 0 1 0 0 0 0 0 90
91 105.5 1 0 0 0 0 0 0 1 0 0 0 0 91
92 98.3 1 0 0 0 0 0 0 0 1 0 0 0 92
93 116.4 1 0 0 0 0 0 0 0 0 1 0 0 93
94 97.4 1 0 0 0 0 0 0 0 0 0 1 0 94
95 93.3 1 0 0 0 0 0 0 0 0 0 0 1 95
96 117.4 1 0 0 0 0 0 0 0 0 0 0 0 96
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Xt_dummy M1 M2 M3 M4
103.8144 -10.6916 3.3810 -5.0478 -13.0892 -10.3441
M5 M6 M7 M8 M9 M10
-9.9355 1.8482 -5.6807 -7.5596 1.9991 -18.5173
M11 t
-20.8961 0.3664
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-20.969 -3.398 0.313 3.964 21.515
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 103.81436 3.42986 30.268 < 2e-16 ***
Xt_dummy -10.69157 2.89231 -3.697 0.000394 ***
M1 3.38103 4.07068 0.831 0.408621
M2 -5.04783 4.06862 -1.241 0.218265
M3 -13.08919 4.06702 -3.218 0.001847 **
M4 -10.34411 4.07175 -2.540 0.012957 *
M5 -9.93547 4.06832 -2.442 0.016750 *
M6 1.84817 4.06535 0.455 0.650587
M7 -5.68069 4.06283 -1.398 0.165821
M8 -7.55955 4.06078 -1.862 0.066242 .
M9 1.99909 4.05917 0.492 0.623692
M10 -18.51728 4.05803 -4.563 1.75e-05 ***
M11 -20.89614 4.05734 -5.150 1.76e-06 ***
t 0.36636 0.04311 8.498 7.18e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.114 on 82 degrees of freedom
Multiple R-squared: 0.669, Adjusted R-squared: 0.6165
F-statistic: 12.75 on 13 and 82 DF, p-value: 9.424e-15
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 9.588258e-03 1.917652e-02 0.9904117
[2,] 1.512951e-03 3.025903e-03 0.9984870
[3,] 2.110355e-04 4.220710e-04 0.9997890
[4,] 4.686202e-03 9.372404e-03 0.9953138
[5,] 2.018180e-03 4.036360e-03 0.9979818
[6,] 5.829844e-04 1.165969e-03 0.9994170
[7,] 1.693182e-04 3.386365e-04 0.9998307
[8,] 4.509214e-05 9.018428e-05 0.9999549
[9,] 1.192426e-05 2.384851e-05 0.9999881
[10,] 3.012720e-06 6.025440e-06 0.9999970
[11,] 4.023824e-05 8.047648e-05 0.9999598
[12,] 1.329317e-05 2.658634e-05 0.9999867
[13,] 4.616250e-06 9.232500e-06 0.9999954
[14,] 7.262562e-05 1.452512e-04 0.9999274
[15,] 3.921933e-05 7.843866e-05 0.9999608
[16,] 1.753969e-05 3.507938e-05 0.9999825
[17,] 1.019436e-04 2.038872e-04 0.9998981
[18,] 4.376016e-05 8.752033e-05 0.9999562
[19,] 3.175774e-05 6.351548e-05 0.9999682
[20,] 2.468629e-05 4.937258e-05 0.9999753
[21,] 1.155879e-05 2.311757e-05 0.9999884
[22,] 7.063087e-06 1.412617e-05 0.9999929
[23,] 7.764594e-06 1.552919e-05 0.9999922
[24,] 5.615171e-06 1.123034e-05 0.9999944
[25,] 3.539623e-06 7.079246e-06 0.9999965
[26,] 3.705986e-06 7.411972e-06 0.9999963
[27,] 1.954604e-06 3.909209e-06 0.9999980
[28,] 9.833851e-07 1.966770e-06 0.9999990
[29,] 7.835294e-07 1.567059e-06 0.9999992
[30,] 1.816839e-06 3.633679e-06 0.9999982
[31,] 1.129149e-06 2.258298e-06 0.9999989
[32,] 6.955803e-07 1.391161e-06 0.9999993
[33,] 1.463663e-05 2.927327e-05 0.9999854
[34,] 1.771523e-05 3.543046e-05 0.9999823
[35,] 1.907436e-05 3.814873e-05 0.9999809
[36,] 1.838892e-05 3.677784e-05 0.9999816
[37,] 1.792212e-05 3.584425e-05 0.9999821
[38,] 3.086057e-05 6.172114e-05 0.9999691
[39,] 6.149311e-05 1.229862e-04 0.9999385
[40,] 1.797417e-04 3.594833e-04 0.9998203
[41,] 2.354224e-04 4.708447e-04 0.9997646
[42,] 1.162713e-03 2.325426e-03 0.9988373
[43,] 2.783979e-03 5.567958e-03 0.9972160
[44,] 1.109757e-02 2.219515e-02 0.9889024
[45,] 1.029322e-01 2.058645e-01 0.8970678
[46,] 1.909115e-01 3.818231e-01 0.8090885
[47,] 4.645295e-01 9.290590e-01 0.5354705
[48,] 4.077844e-01 8.155688e-01 0.5922156
[49,] 3.838581e-01 7.677163e-01 0.6161419
[50,] 3.183776e-01 6.367553e-01 0.6816224
[51,] 4.464692e-01 8.929384e-01 0.5535308
[52,] 4.004639e-01 8.009279e-01 0.5995361
[53,] 3.320571e-01 6.641141e-01 0.6679429
[54,] 3.428736e-01 6.857473e-01 0.6571264
[55,] 2.767330e-01 5.534660e-01 0.7232670
[56,] 5.959215e-01 8.081570e-01 0.4040785
[57,] 6.100915e-01 7.798171e-01 0.3899085
[58,] 6.307602e-01 7.384796e-01 0.3692398
[59,] 5.324511e-01 9.350978e-01 0.4675489
[60,] 4.820055e-01 9.640110e-01 0.5179945
[61,] 3.593000e-01 7.186000e-01 0.6407000
[62,] 4.174547e-01 8.349094e-01 0.5825453
[63,] 5.611872e-01 8.776255e-01 0.4388128
> postscript(file="/var/www/html/rcomp/tmp/1r6vk1261853375.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/2510y1261853375.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/3t2601261853375.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/46s081261853375.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/54s2l1261853375.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 = 96
Frequency = 1
1 2 3 4 5
4.038249196 5.100749196 -0.224250804 3.364303271 1.989303271
6 7 8 9 10
-1.560696729 1.601803271 7.414303271 -1.010696729 4.839303271
11 12 13 14 15
1.251803271 0.689303271 2.241907733 -1.395592267 -2.020592267
16 17 18 19 20
-1.932038192 -1.607038192 -5.757038192 -2.194538192 -6.482038192
21 22 23 24 25
-9.807038192 0.042961808 -5.544538192 -4.007038192 -6.154433730
26 27 28 29 30
-7.091933730 0.583066270 -6.828379656 -5.103379656 2.046620344
31 32 33 34 35
-1.490879656 -6.778379656 1.996620344 -2.953379656 -0.440879656
36 37 38 39 40
0.996620344 -4.850775194 -1.288275194 3.186724806 0.675278881
41 42 43 44 45
0.600278881 3.750278881 1.412778881 -2.874721119 1.900278881
46 47 48 49 50
-9.249721119 -1.737221119 -0.499721119 -13.447116657 -0.984616657
51 52 53 54 55
-0.009616657 -0.821062583 1.003937417 8.253937417 -3.983562583
56 57 58 59 60
9.728937417 4.303937417 -2.546062583 4.366437417 2.703937417
61 62 63 64 65
2.856541879 6.519041879 2.594041879 3.982595954 0.907595954
66 67 68 69 70
3.257595954 -1.579904046 6.432595954 4.307595954 3.957595954
71 72 73 74 75
5.370095954 -1.292404046 4.960200416 3.522700416 -3.202299584
76 77 78 79 80
17.377821894 15.602821894 5.252821894 21.515321894 13.527821894
81 82 83 84 85
11.102821894 17.552821894 10.465321894 12.302821894 10.355426357
86 87 88 89 90
-4.382073643 -0.907073643 -15.818519569 -13.393519569 -15.243519569
91 92 93 94 95
-15.281019569 -20.968519569 -12.793519569 -11.643519569 -13.731019569
96
-10.893519569
> postscript(file="/var/www/html/rcomp/tmp/6rb041261853375.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 = 96
Frequency = 1
lag(myerror, k = 1) myerror
0 4.038249196 NA
1 5.100749196 4.038249196
2 -0.224250804 5.100749196
3 3.364303271 -0.224250804
4 1.989303271 3.364303271
5 -1.560696729 1.989303271
6 1.601803271 -1.560696729
7 7.414303271 1.601803271
8 -1.010696729 7.414303271
9 4.839303271 -1.010696729
10 1.251803271 4.839303271
11 0.689303271 1.251803271
12 2.241907733 0.689303271
13 -1.395592267 2.241907733
14 -2.020592267 -1.395592267
15 -1.932038192 -2.020592267
16 -1.607038192 -1.932038192
17 -5.757038192 -1.607038192
18 -2.194538192 -5.757038192
19 -6.482038192 -2.194538192
20 -9.807038192 -6.482038192
21 0.042961808 -9.807038192
22 -5.544538192 0.042961808
23 -4.007038192 -5.544538192
24 -6.154433730 -4.007038192
25 -7.091933730 -6.154433730
26 0.583066270 -7.091933730
27 -6.828379656 0.583066270
28 -5.103379656 -6.828379656
29 2.046620344 -5.103379656
30 -1.490879656 2.046620344
31 -6.778379656 -1.490879656
32 1.996620344 -6.778379656
33 -2.953379656 1.996620344
34 -0.440879656 -2.953379656
35 0.996620344 -0.440879656
36 -4.850775194 0.996620344
37 -1.288275194 -4.850775194
38 3.186724806 -1.288275194
39 0.675278881 3.186724806
40 0.600278881 0.675278881
41 3.750278881 0.600278881
42 1.412778881 3.750278881
43 -2.874721119 1.412778881
44 1.900278881 -2.874721119
45 -9.249721119 1.900278881
46 -1.737221119 -9.249721119
47 -0.499721119 -1.737221119
48 -13.447116657 -0.499721119
49 -0.984616657 -13.447116657
50 -0.009616657 -0.984616657
51 -0.821062583 -0.009616657
52 1.003937417 -0.821062583
53 8.253937417 1.003937417
54 -3.983562583 8.253937417
55 9.728937417 -3.983562583
56 4.303937417 9.728937417
57 -2.546062583 4.303937417
58 4.366437417 -2.546062583
59 2.703937417 4.366437417
60 2.856541879 2.703937417
61 6.519041879 2.856541879
62 2.594041879 6.519041879
63 3.982595954 2.594041879
64 0.907595954 3.982595954
65 3.257595954 0.907595954
66 -1.579904046 3.257595954
67 6.432595954 -1.579904046
68 4.307595954 6.432595954
69 3.957595954 4.307595954
70 5.370095954 3.957595954
71 -1.292404046 5.370095954
72 4.960200416 -1.292404046
73 3.522700416 4.960200416
74 -3.202299584 3.522700416
75 17.377821894 -3.202299584
76 15.602821894 17.377821894
77 5.252821894 15.602821894
78 21.515321894 5.252821894
79 13.527821894 21.515321894
80 11.102821894 13.527821894
81 17.552821894 11.102821894
82 10.465321894 17.552821894
83 12.302821894 10.465321894
84 10.355426357 12.302821894
85 -4.382073643 10.355426357
86 -0.907073643 -4.382073643
87 -15.818519569 -0.907073643
88 -13.393519569 -15.818519569
89 -15.243519569 -13.393519569
90 -15.281019569 -15.243519569
91 -20.968519569 -15.281019569
92 -12.793519569 -20.968519569
93 -11.643519569 -12.793519569
94 -13.731019569 -11.643519569
95 -10.893519569 -13.731019569
96 NA -10.893519569
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 5.100749196 4.038249196
[2,] -0.224250804 5.100749196
[3,] 3.364303271 -0.224250804
[4,] 1.989303271 3.364303271
[5,] -1.560696729 1.989303271
[6,] 1.601803271 -1.560696729
[7,] 7.414303271 1.601803271
[8,] -1.010696729 7.414303271
[9,] 4.839303271 -1.010696729
[10,] 1.251803271 4.839303271
[11,] 0.689303271 1.251803271
[12,] 2.241907733 0.689303271
[13,] -1.395592267 2.241907733
[14,] -2.020592267 -1.395592267
[15,] -1.932038192 -2.020592267
[16,] -1.607038192 -1.932038192
[17,] -5.757038192 -1.607038192
[18,] -2.194538192 -5.757038192
[19,] -6.482038192 -2.194538192
[20,] -9.807038192 -6.482038192
[21,] 0.042961808 -9.807038192
[22,] -5.544538192 0.042961808
[23,] -4.007038192 -5.544538192
[24,] -6.154433730 -4.007038192
[25,] -7.091933730 -6.154433730
[26,] 0.583066270 -7.091933730
[27,] -6.828379656 0.583066270
[28,] -5.103379656 -6.828379656
[29,] 2.046620344 -5.103379656
[30,] -1.490879656 2.046620344
[31,] -6.778379656 -1.490879656
[32,] 1.996620344 -6.778379656
[33,] -2.953379656 1.996620344
[34,] -0.440879656 -2.953379656
[35,] 0.996620344 -0.440879656
[36,] -4.850775194 0.996620344
[37,] -1.288275194 -4.850775194
[38,] 3.186724806 -1.288275194
[39,] 0.675278881 3.186724806
[40,] 0.600278881 0.675278881
[41,] 3.750278881 0.600278881
[42,] 1.412778881 3.750278881
[43,] -2.874721119 1.412778881
[44,] 1.900278881 -2.874721119
[45,] -9.249721119 1.900278881
[46,] -1.737221119 -9.249721119
[47,] -0.499721119 -1.737221119
[48,] -13.447116657 -0.499721119
[49,] -0.984616657 -13.447116657
[50,] -0.009616657 -0.984616657
[51,] -0.821062583 -0.009616657
[52,] 1.003937417 -0.821062583
[53,] 8.253937417 1.003937417
[54,] -3.983562583 8.253937417
[55,] 9.728937417 -3.983562583
[56,] 4.303937417 9.728937417
[57,] -2.546062583 4.303937417
[58,] 4.366437417 -2.546062583
[59,] 2.703937417 4.366437417
[60,] 2.856541879 2.703937417
[61,] 6.519041879 2.856541879
[62,] 2.594041879 6.519041879
[63,] 3.982595954 2.594041879
[64,] 0.907595954 3.982595954
[65,] 3.257595954 0.907595954
[66,] -1.579904046 3.257595954
[67,] 6.432595954 -1.579904046
[68,] 4.307595954 6.432595954
[69,] 3.957595954 4.307595954
[70,] 5.370095954 3.957595954
[71,] -1.292404046 5.370095954
[72,] 4.960200416 -1.292404046
[73,] 3.522700416 4.960200416
[74,] -3.202299584 3.522700416
[75,] 17.377821894 -3.202299584
[76,] 15.602821894 17.377821894
[77,] 5.252821894 15.602821894
[78,] 21.515321894 5.252821894
[79,] 13.527821894 21.515321894
[80,] 11.102821894 13.527821894
[81,] 17.552821894 11.102821894
[82,] 10.465321894 17.552821894
[83,] 12.302821894 10.465321894
[84,] 10.355426357 12.302821894
[85,] -4.382073643 10.355426357
[86,] -0.907073643 -4.382073643
[87,] -15.818519569 -0.907073643
[88,] -13.393519569 -15.818519569
[89,] -15.243519569 -13.393519569
[90,] -15.281019569 -15.243519569
[91,] -20.968519569 -15.281019569
[92,] -12.793519569 -20.968519569
[93,] -11.643519569 -12.793519569
[94,] -13.731019569 -11.643519569
[95,] -10.893519569 -13.731019569
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 5.100749196 4.038249196
2 -0.224250804 5.100749196
3 3.364303271 -0.224250804
4 1.989303271 3.364303271
5 -1.560696729 1.989303271
6 1.601803271 -1.560696729
7 7.414303271 1.601803271
8 -1.010696729 7.414303271
9 4.839303271 -1.010696729
10 1.251803271 4.839303271
11 0.689303271 1.251803271
12 2.241907733 0.689303271
13 -1.395592267 2.241907733
14 -2.020592267 -1.395592267
15 -1.932038192 -2.020592267
16 -1.607038192 -1.932038192
17 -5.757038192 -1.607038192
18 -2.194538192 -5.757038192
19 -6.482038192 -2.194538192
20 -9.807038192 -6.482038192
21 0.042961808 -9.807038192
22 -5.544538192 0.042961808
23 -4.007038192 -5.544538192
24 -6.154433730 -4.007038192
25 -7.091933730 -6.154433730
26 0.583066270 -7.091933730
27 -6.828379656 0.583066270
28 -5.103379656 -6.828379656
29 2.046620344 -5.103379656
30 -1.490879656 2.046620344
31 -6.778379656 -1.490879656
32 1.996620344 -6.778379656
33 -2.953379656 1.996620344
34 -0.440879656 -2.953379656
35 0.996620344 -0.440879656
36 -4.850775194 0.996620344
37 -1.288275194 -4.850775194
38 3.186724806 -1.288275194
39 0.675278881 3.186724806
40 0.600278881 0.675278881
41 3.750278881 0.600278881
42 1.412778881 3.750278881
43 -2.874721119 1.412778881
44 1.900278881 -2.874721119
45 -9.249721119 1.900278881
46 -1.737221119 -9.249721119
47 -0.499721119 -1.737221119
48 -13.447116657 -0.499721119
49 -0.984616657 -13.447116657
50 -0.009616657 -0.984616657
51 -0.821062583 -0.009616657
52 1.003937417 -0.821062583
53 8.253937417 1.003937417
54 -3.983562583 8.253937417
55 9.728937417 -3.983562583
56 4.303937417 9.728937417
57 -2.546062583 4.303937417
58 4.366437417 -2.546062583
59 2.703937417 4.366437417
60 2.856541879 2.703937417
61 6.519041879 2.856541879
62 2.594041879 6.519041879
63 3.982595954 2.594041879
64 0.907595954 3.982595954
65 3.257595954 0.907595954
66 -1.579904046 3.257595954
67 6.432595954 -1.579904046
68 4.307595954 6.432595954
69 3.957595954 4.307595954
70 5.370095954 3.957595954
71 -1.292404046 5.370095954
72 4.960200416 -1.292404046
73 3.522700416 4.960200416
74 -3.202299584 3.522700416
75 17.377821894 -3.202299584
76 15.602821894 17.377821894
77 5.252821894 15.602821894
78 21.515321894 5.252821894
79 13.527821894 21.515321894
80 11.102821894 13.527821894
81 17.552821894 11.102821894
82 10.465321894 17.552821894
83 12.302821894 10.465321894
84 10.355426357 12.302821894
85 -4.382073643 10.355426357
86 -0.907073643 -4.382073643
87 -15.818519569 -0.907073643
88 -13.393519569 -15.818519569
89 -15.243519569 -13.393519569
90 -15.281019569 -15.243519569
91 -20.968519569 -15.281019569
92 -12.793519569 -20.968519569
93 -11.643519569 -12.793519569
94 -13.731019569 -11.643519569
95 -10.893519569 -13.731019569
> 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/76h691261853375.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/8c8lu1261853375.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/9gyex1261853375.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/101stp1261853375.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/1192981261853375.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/12vg871261853375.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/131xqc1261853375.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/14sv7y1261853375.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/151ziq1261853375.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/16qh6j1261853375.tab")
+ }
>
> try(system("convert tmp/1r6vk1261853375.ps tmp/1r6vk1261853375.png",intern=TRUE))
character(0)
> try(system("convert tmp/2510y1261853375.ps tmp/2510y1261853375.png",intern=TRUE))
character(0)
> try(system("convert tmp/3t2601261853375.ps tmp/3t2601261853375.png",intern=TRUE))
character(0)
> try(system("convert tmp/46s081261853375.ps tmp/46s081261853375.png",intern=TRUE))
character(0)
> try(system("convert tmp/54s2l1261853375.ps tmp/54s2l1261853375.png",intern=TRUE))
character(0)
> try(system("convert tmp/6rb041261853375.ps tmp/6rb041261853375.png",intern=TRUE))
character(0)
> try(system("convert tmp/76h691261853375.ps tmp/76h691261853375.png",intern=TRUE))
character(0)
> try(system("convert tmp/8c8lu1261853375.ps tmp/8c8lu1261853375.png",intern=TRUE))
character(0)
> try(system("convert tmp/9gyex1261853375.ps tmp/9gyex1261853375.png",intern=TRUE))
character(0)
> try(system("convert tmp/101stp1261853375.ps tmp/101stp1261853375.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.961 1.610 4.887