R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-pc-linux-gnu (32-bit)
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(104.29
+ ,103.65
+ ,104.12
+ ,106.67
+ ,105.03
+ ,104.56
+ ,103.87
+ ,104.76
+ ,106.86
+ ,105.32
+ ,104.79
+ ,103.94
+ ,105.37
+ ,107.22
+ ,105.52
+ ,105.08
+ ,105.32
+ ,104.97
+ ,107.5
+ ,105.67
+ ,105.21
+ ,105.54
+ ,105.63
+ ,107.35
+ ,105.71
+ ,105.43
+ ,106.08
+ ,106.17
+ ,107.45
+ ,105.81
+ ,105.69
+ ,106.21
+ ,106.05
+ ,108.23
+ ,105.99
+ ,105.74
+ ,105.53
+ ,106.21
+ ,108.39
+ ,106.02
+ ,106.2
+ ,105.56
+ ,108.06
+ ,108
+ ,106.19
+ ,106.04
+ ,105.14
+ ,107.95
+ ,107.59
+ ,106.22
+ ,106.45
+ ,105.97
+ ,108.22
+ ,107.74
+ ,106.34
+ ,106.4
+ ,105.45
+ ,107.56
+ ,108.17
+ ,106.42
+ ,106.48
+ ,106.22
+ ,106.7
+ ,108.44
+ ,106.84
+ ,106.83
+ ,106.31
+ ,107.38
+ ,108.85
+ ,107.23
+ ,107.14
+ ,107.37
+ ,107.42
+ ,108.8
+ ,107.42
+ ,107.94
+ ,109.31
+ ,108.17
+ ,109.46
+ ,107.63
+ ,108.46
+ ,110.82
+ ,108.89
+ ,109.56
+ ,107.69
+ ,108.81
+ ,111.22
+ ,108.87
+ ,109.94
+ ,107.81
+ ,108.92
+ ,110.66
+ ,108.24
+ ,111.06
+ ,107.92
+ ,108.99
+ ,110.76
+ ,108.23
+ ,110.9
+ ,108.06
+ ,109.16
+ ,110.69
+ ,109.03
+ ,110.79
+ ,108.21
+ ,109.22
+ ,111.08
+ ,108.24
+ ,111.08
+ ,108.44
+ ,109.43
+ ,110.97
+ ,108.01
+ ,111.91
+ ,108.55
+ ,109.23
+ ,110.24
+ ,107.72
+ ,112.09
+ ,108.66
+ ,109.93
+ ,112.51
+ ,107.81
+ ,112.43
+ ,109.23
+ ,110.09
+ ,111.52
+ ,107.98
+ ,113.44
+ ,109.7
+ ,110.33
+ ,112.13
+ ,108.34
+ ,113.4
+ ,109.94
+ ,110.11
+ ,112.23
+ ,108.91
+ ,112.5
+ ,110.13
+ ,110.35
+ ,112.92
+ ,108.78
+ ,112.73
+ ,110.39
+ ,110.09
+ ,111.89
+ ,108.34
+ ,113.12
+ ,110.46
+ ,110.44
+ ,111.99
+ ,108.64
+ ,113.77
+ ,110.67
+ ,110.39
+ ,111.51
+ ,108.68
+ ,113.93
+ ,110.89
+ ,110.62
+ ,112.33
+ ,109.31
+ ,113.41
+ ,110.98
+ ,110.43
+ ,112.04
+ ,109.65
+ ,112.62
+ ,111.12
+ ,110.46
+ ,112.09
+ ,109.07
+ ,113.12
+ ,111.33
+ ,110.55
+ ,111.41
+ ,109.18
+ ,113.65
+ ,111.43
+ ,110.94
+ ,112.61
+ ,109.71
+ ,113.55
+ ,111.87
+ ,111.56
+ ,113.14
+ ,110.68
+ ,114.28
+ ,112.22
+ ,111.82
+ ,113.65
+ ,111.09
+ ,114.31
+ ,112.47
+ ,111.73
+ ,114.26
+ ,109.64
+ ,115.09
+ ,112.64
+ ,111.57
+ ,114.4
+ ,109.08
+ ,114.73
+ ,112.84
+ ,111.85
+ ,114.93
+ ,109.27
+ ,115.13
+ ,113.03
+ ,112.06
+ ,114.86
+ ,109.41
+ ,115.74
+ ,113.09
+ ,112.2
+ ,114.95
+ ,109.99
+ ,115.78
+ ,113.27
+ ,112.47
+ ,116.17
+ ,110.35
+ ,115.42
+ ,113.44
+ ,112.15
+ ,114.6
+ ,110.25
+ ,115.44
+ ,113.51
+ ,112.36
+ ,114.62
+ ,110.33
+ ,116
+ ,113.66
+ ,112.32
+ ,113.82
+ ,110.29
+ ,116.44
+ ,113.62
+ ,112.67
+ ,115.02
+ ,110.45
+ ,116.38
+ ,114.01
+ ,113.02
+ ,115.18
+ ,110.75
+ ,117.17
+ ,114.55
+ ,113.05
+ ,115.59
+ ,111.15
+ ,116.75
+ ,114.77
+ ,113.5
+ ,116.6
+ ,111.56
+ ,117.5
+ ,114.87
+ ,113.67
+ ,117.07
+ ,112.33
+ ,117.43
+ ,115.11
+ ,113.65
+ ,116.96
+ ,112.13
+ ,117.65
+ ,115.09
+ ,114
+ ,116.66
+ ,112.49
+ ,118.65
+ ,115.24
+ ,114.03
+ ,116.07
+ ,113.14
+ ,118.58
+ ,115.27
+ ,114.08
+ ,116.04
+ ,113.42
+ ,118.42
+ ,115.41
+ ,114.49
+ ,115.81
+ ,114.67
+ ,118.55
+ ,115.59
+ ,114.48
+ ,116.22
+ ,114.03
+ ,118.77
+ ,115.6
+ ,114.25
+ ,115.85
+ ,113.37
+ ,118.71
+ ,115.68
+ ,114.68
+ ,116.43
+ ,113.2
+ ,119.58
+ ,116.19
+ ,115.28
+ ,117.39
+ ,114.2
+ ,119.97
+ ,116.55
+ ,115.9
+ ,119.17
+ ,114.97
+ ,119.99
+ ,116.73
+ ,115.87
+ ,119.24
+ ,115.72
+ ,119.67
+ ,117.04
+ ,116.09
+ ,120.03
+ ,115.47
+ ,120.04
+ ,117.12
+ ,116.29
+ ,119.34
+ ,116.3
+ ,120.51
+ ,117.28
+ ,116.76
+ ,118.49
+ ,117.66
+ ,121.47
+ ,117.48
+ ,116.78
+ ,118.59
+ ,118.01
+ ,121.2
+ ,117.66
+ ,116.65
+ ,117.5
+ ,119.07
+ ,120.81
+ ,117.92
+ ,116.46
+ ,117.56
+ ,118.29
+ ,121.19
+ ,118.12
+ ,116.82
+ ,118.25
+ ,117.57
+ ,121.67
+ ,118.17
+ ,116.91
+ ,118.01
+ ,117.61
+ ,121.67
+ ,118.39)
+ ,dim=c(5
+ ,72)
+ ,dimnames=list(c('index'
+ ,'voeding'
+ ,'nietvoeding'
+ ,'diensten'
+ ,'huur')
+ ,1:72))
> y <- array(NA,dim=c(5,72),dimnames=list(c('index','voeding','nietvoeding','diensten','huur'),1:72))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Do not include Seasonal Dummies'
> par1 = '1'
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from 'package:base':
as.Date, 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
index voeding nietvoeding diensten huur
1 104.29 103.65 104.12 106.67 105.03
2 104.56 103.87 104.76 106.86 105.32
3 104.79 103.94 105.37 107.22 105.52
4 105.08 105.32 104.97 107.50 105.67
5 105.21 105.54 105.63 107.35 105.71
6 105.43 106.08 106.17 107.45 105.81
7 105.69 106.21 106.05 108.23 105.99
8 105.74 105.53 106.21 108.39 106.02
9 106.20 105.56 108.06 108.00 106.19
10 106.04 105.14 107.95 107.59 106.22
11 106.45 105.97 108.22 107.74 106.34
12 106.40 105.45 107.56 108.17 106.42
13 106.48 106.22 106.70 108.44 106.84
14 106.83 106.31 107.38 108.85 107.23
15 107.14 107.37 107.42 108.80 107.42
16 107.94 109.31 108.17 109.46 107.63
17 108.46 110.82 108.89 109.56 107.69
18 108.81 111.22 108.87 109.94 107.81
19 108.92 110.66 108.24 111.06 107.92
20 108.99 110.76 108.23 110.90 108.06
21 109.16 110.69 109.03 110.79 108.21
22 109.22 111.08 108.24 111.08 108.44
23 109.43 110.97 108.01 111.91 108.55
24 109.23 110.24 107.72 112.09 108.66
25 109.93 112.51 107.81 112.43 109.23
26 110.09 111.52 107.98 113.44 109.70
27 110.33 112.13 108.34 113.40 109.94
28 110.11 112.23 108.91 112.50 110.13
29 110.35 112.92 108.78 112.73 110.39
30 110.09 111.89 108.34 113.12 110.46
31 110.44 111.99 108.64 113.77 110.67
32 110.39 111.51 108.68 113.93 110.89
33 110.62 112.33 109.31 113.41 110.98
34 110.43 112.04 109.65 112.62 111.12
35 110.46 112.09 109.07 113.12 111.33
36 110.55 111.41 109.18 113.65 111.43
37 110.94 112.61 109.71 113.55 111.87
38 111.56 113.14 110.68 114.28 112.22
39 111.82 113.65 111.09 114.31 112.47
40 111.73 114.26 109.64 115.09 112.64
41 111.57 114.40 109.08 114.73 112.84
42 111.85 114.93 109.27 115.13 113.03
43 112.06 114.86 109.41 115.74 113.09
44 112.20 114.95 109.99 115.78 113.27
45 112.47 116.17 110.35 115.42 113.44
46 112.15 114.60 110.25 115.44 113.51
47 112.36 114.62 110.33 116.00 113.66
48 112.32 113.82 110.29 116.44 113.62
49 112.67 115.02 110.45 116.38 114.01
50 113.02 115.18 110.75 117.17 114.55
51 113.05 115.59 111.15 116.75 114.77
52 113.50 116.60 111.56 117.50 114.87
53 113.67 117.07 112.33 117.43 115.11
54 113.65 116.96 112.13 117.65 115.09
55 114.00 116.66 112.49 118.65 115.24
56 114.03 116.07 113.14 118.58 115.27
57 114.08 116.04 113.42 118.42 115.41
58 114.49 115.81 114.67 118.55 115.59
59 114.48 116.22 114.03 118.77 115.60
60 114.25 115.85 113.37 118.71 115.68
61 114.68 116.43 113.20 119.58 116.19
62 115.28 117.39 114.20 119.97 116.55
63 115.90 119.17 114.97 119.99 116.73
64 115.87 119.24 115.72 119.67 117.04
65 116.09 120.03 115.47 120.04 117.12
66 116.29 119.34 116.30 120.51 117.28
67 116.76 118.49 117.66 121.47 117.48
68 116.78 118.59 118.01 121.20 117.66
69 116.65 117.50 119.07 120.81 117.92
70 116.46 117.56 118.29 121.19 118.12
71 116.82 118.25 117.57 121.67 118.17
72 116.91 118.01 117.61 121.67 118.39
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) voeding nietvoeding diensten huur
14.76361 0.30686 0.19074 0.33460 0.02206
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.44550 -0.10084 0.02433 0.12597 0.28826
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.76361 0.76224 19.369 < 2e-16 ***
voeding 0.30686 0.01851 16.580 < 2e-16 ***
nietvoeding 0.19074 0.01610 11.850 < 2e-16 ***
diensten 0.33460 0.04094 8.172 1.18e-11 ***
huur 0.02206 0.04342 0.508 0.613
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1748 on 67 degrees of freedom
Multiple R-squared: 0.9977, Adjusted R-squared: 0.9976
F-statistic: 7402 on 4 and 67 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.03985832 0.079716634 0.960141683
[2,] 0.12201414 0.244028274 0.877985863
[3,] 0.09862010 0.197240196 0.901379902
[4,] 0.30458076 0.609161529 0.695419236
[5,] 0.38492304 0.769846089 0.615076955
[6,] 0.34870720 0.697414397 0.651292801
[7,] 0.44541389 0.890827783 0.554586108
[8,] 0.43283099 0.865661972 0.567169014
[9,] 0.68479006 0.630419880 0.315209940
[10,] 0.83154477 0.336910456 0.168455228
[11,] 0.93978911 0.120421780 0.060210890
[12,] 0.97480193 0.050396144 0.025198072
[13,] 0.97586236 0.048275285 0.024137643
[14,] 0.96669110 0.066617794 0.033308897
[15,] 0.95053500 0.098929995 0.049464998
[16,] 0.92872777 0.142544458 0.071272229
[17,] 0.91431166 0.171376673 0.085688336
[18,] 0.91989793 0.160204148 0.080102074
[19,] 0.96710997 0.065780058 0.032890029
[20,] 0.98362798 0.032744047 0.016372023
[21,] 0.99517555 0.009648896 0.004824448
[22,] 0.99560282 0.008794361 0.004397180
[23,] 0.99645660 0.007086800 0.003543400
[24,] 0.99635233 0.007295341 0.003647670
[25,] 0.99636843 0.007263141 0.003631570
[26,] 0.99451023 0.010979542 0.005489771
[27,] 0.99106577 0.017868463 0.008934231
[28,] 0.98626260 0.027474798 0.013737399
[29,] 0.97849012 0.043019767 0.021509883
[30,] 0.96739650 0.065206990 0.032603495
[31,] 0.96279370 0.074412591 0.037206295
[32,] 0.97207189 0.055856226 0.027928113
[33,] 0.97356648 0.052867046 0.026433523
[34,] 0.96269146 0.074617089 0.037308544
[35,] 0.95231149 0.095377015 0.047688508
[36,] 0.94726877 0.105462451 0.052731225
[37,] 0.94866938 0.102661243 0.051330621
[38,] 0.94886373 0.102272537 0.051136268
[39,] 0.93884776 0.122304475 0.061152237
[40,] 0.93441674 0.131166522 0.065583261
[41,] 0.93043355 0.139132894 0.069566447
[42,] 0.95313857 0.093722856 0.046861428
[43,] 0.95624986 0.087500284 0.043750142
[44,] 0.96682375 0.066352492 0.033176246
[45,] 0.97544017 0.049119664 0.024559832
[46,] 0.98497599 0.030048012 0.015024006
[47,] 0.99103004 0.017939924 0.008969962
[48,] 0.99682664 0.006346716 0.003173358
[49,] 0.99782665 0.004346706 0.002173353
[50,] 0.99745999 0.005080022 0.002540011
[51,] 0.99521220 0.009575609 0.004787805
[52,] 0.98854582 0.022908358 0.011454179
[53,] 0.97377105 0.052457900 0.026228950
[54,] 0.94597651 0.108046974 0.054023487
[55,] 0.89952958 0.200940839 0.100470420
[56,] 0.85060723 0.298785534 0.149392767
[57,] 0.81697527 0.366049459 0.183024729
> postscript(file="/var/fisher/rcomp/tmp/1qhbn1352991132.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/28z581352991132.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/3yvat1352991132.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/4f4ny1352991132.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/5v3ao1352991132.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> qqnorm(mysum$resid, main='Residual Normal Q-Q Plot')
> qqline(mysum$resid)
> grid()
> dev.off()
null device
1
> (myerror <- as.ts(mysum$resid))
Time Series:
Start = 1
End = 72
Frequency = 1
1 2 3 4 5 6
-0.148649392 -0.138205367 -0.170907691 -0.325073877 -0.339164846 -0.423536071
7 8 9 10 11 12
-0.445500231 -0.271553324 -0.046887597 0.079500872 0.130469728 0.220282319
13 14 15 16 17 18
0.128432289 0.175319587 0.164958490 0.001123921 -0.114351986 -0.013077258
19 20 21 22 23 24
0.011748834 0.103418917 0.175802546 0.164705867 0.172183438 0.188850870
25 26 27 28 29 30
0.048775585 0.131822417 0.124061249 0.061604633 0.031974468 0.039926493
31 32 33 34 35 36
0.079892908 0.111166154 0.141382669 0.236768379 0.190122072 0.288258549
37 38 39 40 41 42
0.232688931 0.253051815 0.262796542 -0.002551394 0.017350247 -0.039558832
43 44 45 46 47 48
-0.040214680 -0.055817425 -0.112145250 0.060461619 0.058378027 0.125151487
49 50 51 52 53 54
0.087875451 0.055306656 0.018878306 -0.172412780 -0.275380126 -0.296648809
55 56 57 58 59 60
-0.261171313 -0.151346771 -0.095100437 0.099579635 0.012009362 0.039749112
61 62 63 64 65 66
0.021842060 -0.001921914 -0.085665521 -0.179967408 -0.280268500 -0.187645509
67 68 69 70 71 72
-0.041857111 -0.032930315 0.094119198 -0.097074056 0.026815179 0.177978912
> postscript(file="/var/fisher/rcomp/tmp/63nwf1352991132.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 72
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.148649392 NA
1 -0.138205367 -0.148649392
2 -0.170907691 -0.138205367
3 -0.325073877 -0.170907691
4 -0.339164846 -0.325073877
5 -0.423536071 -0.339164846
6 -0.445500231 -0.423536071
7 -0.271553324 -0.445500231
8 -0.046887597 -0.271553324
9 0.079500872 -0.046887597
10 0.130469728 0.079500872
11 0.220282319 0.130469728
12 0.128432289 0.220282319
13 0.175319587 0.128432289
14 0.164958490 0.175319587
15 0.001123921 0.164958490
16 -0.114351986 0.001123921
17 -0.013077258 -0.114351986
18 0.011748834 -0.013077258
19 0.103418917 0.011748834
20 0.175802546 0.103418917
21 0.164705867 0.175802546
22 0.172183438 0.164705867
23 0.188850870 0.172183438
24 0.048775585 0.188850870
25 0.131822417 0.048775585
26 0.124061249 0.131822417
27 0.061604633 0.124061249
28 0.031974468 0.061604633
29 0.039926493 0.031974468
30 0.079892908 0.039926493
31 0.111166154 0.079892908
32 0.141382669 0.111166154
33 0.236768379 0.141382669
34 0.190122072 0.236768379
35 0.288258549 0.190122072
36 0.232688931 0.288258549
37 0.253051815 0.232688931
38 0.262796542 0.253051815
39 -0.002551394 0.262796542
40 0.017350247 -0.002551394
41 -0.039558832 0.017350247
42 -0.040214680 -0.039558832
43 -0.055817425 -0.040214680
44 -0.112145250 -0.055817425
45 0.060461619 -0.112145250
46 0.058378027 0.060461619
47 0.125151487 0.058378027
48 0.087875451 0.125151487
49 0.055306656 0.087875451
50 0.018878306 0.055306656
51 -0.172412780 0.018878306
52 -0.275380126 -0.172412780
53 -0.296648809 -0.275380126
54 -0.261171313 -0.296648809
55 -0.151346771 -0.261171313
56 -0.095100437 -0.151346771
57 0.099579635 -0.095100437
58 0.012009362 0.099579635
59 0.039749112 0.012009362
60 0.021842060 0.039749112
61 -0.001921914 0.021842060
62 -0.085665521 -0.001921914
63 -0.179967408 -0.085665521
64 -0.280268500 -0.179967408
65 -0.187645509 -0.280268500
66 -0.041857111 -0.187645509
67 -0.032930315 -0.041857111
68 0.094119198 -0.032930315
69 -0.097074056 0.094119198
70 0.026815179 -0.097074056
71 0.177978912 0.026815179
72 NA 0.177978912
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.138205367 -0.148649392
[2,] -0.170907691 -0.138205367
[3,] -0.325073877 -0.170907691
[4,] -0.339164846 -0.325073877
[5,] -0.423536071 -0.339164846
[6,] -0.445500231 -0.423536071
[7,] -0.271553324 -0.445500231
[8,] -0.046887597 -0.271553324
[9,] 0.079500872 -0.046887597
[10,] 0.130469728 0.079500872
[11,] 0.220282319 0.130469728
[12,] 0.128432289 0.220282319
[13,] 0.175319587 0.128432289
[14,] 0.164958490 0.175319587
[15,] 0.001123921 0.164958490
[16,] -0.114351986 0.001123921
[17,] -0.013077258 -0.114351986
[18,] 0.011748834 -0.013077258
[19,] 0.103418917 0.011748834
[20,] 0.175802546 0.103418917
[21,] 0.164705867 0.175802546
[22,] 0.172183438 0.164705867
[23,] 0.188850870 0.172183438
[24,] 0.048775585 0.188850870
[25,] 0.131822417 0.048775585
[26,] 0.124061249 0.131822417
[27,] 0.061604633 0.124061249
[28,] 0.031974468 0.061604633
[29,] 0.039926493 0.031974468
[30,] 0.079892908 0.039926493
[31,] 0.111166154 0.079892908
[32,] 0.141382669 0.111166154
[33,] 0.236768379 0.141382669
[34,] 0.190122072 0.236768379
[35,] 0.288258549 0.190122072
[36,] 0.232688931 0.288258549
[37,] 0.253051815 0.232688931
[38,] 0.262796542 0.253051815
[39,] -0.002551394 0.262796542
[40,] 0.017350247 -0.002551394
[41,] -0.039558832 0.017350247
[42,] -0.040214680 -0.039558832
[43,] -0.055817425 -0.040214680
[44,] -0.112145250 -0.055817425
[45,] 0.060461619 -0.112145250
[46,] 0.058378027 0.060461619
[47,] 0.125151487 0.058378027
[48,] 0.087875451 0.125151487
[49,] 0.055306656 0.087875451
[50,] 0.018878306 0.055306656
[51,] -0.172412780 0.018878306
[52,] -0.275380126 -0.172412780
[53,] -0.296648809 -0.275380126
[54,] -0.261171313 -0.296648809
[55,] -0.151346771 -0.261171313
[56,] -0.095100437 -0.151346771
[57,] 0.099579635 -0.095100437
[58,] 0.012009362 0.099579635
[59,] 0.039749112 0.012009362
[60,] 0.021842060 0.039749112
[61,] -0.001921914 0.021842060
[62,] -0.085665521 -0.001921914
[63,] -0.179967408 -0.085665521
[64,] -0.280268500 -0.179967408
[65,] -0.187645509 -0.280268500
[66,] -0.041857111 -0.187645509
[67,] -0.032930315 -0.041857111
[68,] 0.094119198 -0.032930315
[69,] -0.097074056 0.094119198
[70,] 0.026815179 -0.097074056
[71,] 0.177978912 0.026815179
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.138205367 -0.148649392
2 -0.170907691 -0.138205367
3 -0.325073877 -0.170907691
4 -0.339164846 -0.325073877
5 -0.423536071 -0.339164846
6 -0.445500231 -0.423536071
7 -0.271553324 -0.445500231
8 -0.046887597 -0.271553324
9 0.079500872 -0.046887597
10 0.130469728 0.079500872
11 0.220282319 0.130469728
12 0.128432289 0.220282319
13 0.175319587 0.128432289
14 0.164958490 0.175319587
15 0.001123921 0.164958490
16 -0.114351986 0.001123921
17 -0.013077258 -0.114351986
18 0.011748834 -0.013077258
19 0.103418917 0.011748834
20 0.175802546 0.103418917
21 0.164705867 0.175802546
22 0.172183438 0.164705867
23 0.188850870 0.172183438
24 0.048775585 0.188850870
25 0.131822417 0.048775585
26 0.124061249 0.131822417
27 0.061604633 0.124061249
28 0.031974468 0.061604633
29 0.039926493 0.031974468
30 0.079892908 0.039926493
31 0.111166154 0.079892908
32 0.141382669 0.111166154
33 0.236768379 0.141382669
34 0.190122072 0.236768379
35 0.288258549 0.190122072
36 0.232688931 0.288258549
37 0.253051815 0.232688931
38 0.262796542 0.253051815
39 -0.002551394 0.262796542
40 0.017350247 -0.002551394
41 -0.039558832 0.017350247
42 -0.040214680 -0.039558832
43 -0.055817425 -0.040214680
44 -0.112145250 -0.055817425
45 0.060461619 -0.112145250
46 0.058378027 0.060461619
47 0.125151487 0.058378027
48 0.087875451 0.125151487
49 0.055306656 0.087875451
50 0.018878306 0.055306656
51 -0.172412780 0.018878306
52 -0.275380126 -0.172412780
53 -0.296648809 -0.275380126
54 -0.261171313 -0.296648809
55 -0.151346771 -0.261171313
56 -0.095100437 -0.151346771
57 0.099579635 -0.095100437
58 0.012009362 0.099579635
59 0.039749112 0.012009362
60 0.021842060 0.039749112
61 -0.001921914 0.021842060
62 -0.085665521 -0.001921914
63 -0.179967408 -0.085665521
64 -0.280268500 -0.179967408
65 -0.187645509 -0.280268500
66 -0.041857111 -0.187645509
67 -0.032930315 -0.041857111
68 0.094119198 -0.032930315
69 -0.097074056 0.094119198
70 0.026815179 -0.097074056
71 0.177978912 0.026815179
> 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/fisher/rcomp/tmp/7hnzv1352991132.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/894q81352991132.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/91run1352991132.ps",horizontal=F,onefile=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/fisher/rcomp/tmp/10mls61352991132.ps",horizontal=F,onefile=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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/fisher/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/fisher/rcomp/tmp/11efcd1352991132.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/fisher/rcomp/tmp/12kfmj1352991132.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/fisher/rcomp/tmp/138ojk1352991132.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/fisher/rcomp/tmp/14av0b1352991133.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/fisher/rcomp/tmp/15tdra1352991133.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/fisher/rcomp/tmp/1699wo1352991133.tab")
+ }
>
> try(system("convert tmp/1qhbn1352991132.ps tmp/1qhbn1352991132.png",intern=TRUE))
character(0)
> try(system("convert tmp/28z581352991132.ps tmp/28z581352991132.png",intern=TRUE))
character(0)
> try(system("convert tmp/3yvat1352991132.ps tmp/3yvat1352991132.png",intern=TRUE))
character(0)
> try(system("convert tmp/4f4ny1352991132.ps tmp/4f4ny1352991132.png",intern=TRUE))
character(0)
> try(system("convert tmp/5v3ao1352991132.ps tmp/5v3ao1352991132.png",intern=TRUE))
character(0)
> try(system("convert tmp/63nwf1352991132.ps tmp/63nwf1352991132.png",intern=TRUE))
character(0)
> try(system("convert tmp/7hnzv1352991132.ps tmp/7hnzv1352991132.png",intern=TRUE))
character(0)
> try(system("convert tmp/894q81352991132.ps tmp/894q81352991132.png",intern=TRUE))
character(0)
> try(system("convert tmp/91run1352991132.ps tmp/91run1352991132.png",intern=TRUE))
character(0)
> try(system("convert tmp/10mls61352991132.ps tmp/10mls61352991132.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.706 1.349 8.050