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(105.3
+ ,97.2
+ ,95
+ ,109.3
+ ,108
+ ,97.4
+ ,98.5
+ ,110
+ ,100.9
+ ,101.4
+ ,115.7
+ ,113.8
+ ,100.9
+ ,98
+ ,129.2
+ ,109.4
+ ,112.8
+ ,118.5
+ ,122.1
+ ,109.2
+ ,99.9
+ ,99.9
+ ,101.3
+ ,76.4
+ ,100.9
+ ,100.6
+ ,103.1
+ ,98.5
+ ,100.6
+ ,101.2
+ ,93.5
+ ,99.9
+ ,100.2
+ ,101.8
+ ,104.7
+ ,117.2
+ ,104.8
+ ,117.6
+ ,109.9
+ ,112.3
+ ,103.9
+ ,100.4
+ ,136
+ ,111.4
+ ,118.6
+ ,120.7
+ ,125.7
+ ,110.9
+ ,108.6
+ ,126
+ ,106.9
+ ,110.1
+ ,113.3
+ ,117.5
+ ,106.6
+ ,109.9
+ ,125
+ ,103.9
+ ,113.1
+ ,96.4
+ ,105.8
+ ,103.2
+ ,109.5
+ ,112.5
+ ,107.6
+ ,91.2
+ ,116.2
+ ,115.1
+ ,108.7
+ ,109.2
+ ,119.3
+ ,104.7
+ ,106.2
+ ,119.2
+ ,119.3
+ ,104.6
+ ,100.9
+ ,121.1
+ ,109.1
+ ,115.5
+ ,117.9
+ ,119
+ ,108.7
+ ,98.4
+ ,127.1
+ ,109.3
+ ,106.2
+ ,126.5
+ ,126.8
+ ,109.5
+ ,94.2
+ ,116
+ ,102.2
+ ,95.9
+ ,116.3
+ ,116.3
+ ,102.6
+ ,94.7
+ ,131.9
+ ,109.8
+ ,113.2
+ ,124.7
+ ,127.1
+ ,109.5
+ ,95.2
+ ,101
+ ,106.2
+ ,78.3
+ ,108.8
+ ,106.5
+ ,108.1
+ ,100.3
+ ,92.6
+ ,95.1
+ ,79.8
+ ,96.7
+ ,95.5
+ ,96.1
+ ,100.9
+ ,127.2
+ ,118.7
+ ,121.2
+ ,118.4
+ ,121.3
+ ,118.5
+ ,97.9
+ ,124.3
+ ,116.9
+ ,125.6
+ ,115.4
+ ,118.3
+ ,116.3
+ ,106.9
+ ,103.8
+ ,105.3
+ ,97.2
+ ,91.7
+ ,95.6
+ ,105.9
+ ,100.8
+ ,106.4
+ ,119.5
+ ,102.8
+ ,82.5
+ ,90.1
+ ,120.7
+ ,106.6
+ ,84.2
+ ,96.5
+ ,88.8
+ ,81.7
+ ,82.6
+ ,97
+ ,108.2
+ ,91.9
+ ,99.3
+ ,95.3
+ ,84.5
+ ,86.9
+ ,99.6
+ ,98.4
+ ,103.4
+ ,113.8
+ ,107.6
+ ,92.4
+ ,95.9
+ ,114.2
+ ,102
+ ,93
+ ,102.7
+ ,95
+ ,86.8
+ ,88.8
+ ,103.3
+ ,95.7
+ ,110.6
+ ,98.8
+ ,87.5
+ ,85.1
+ ,93.2
+ ,99.6
+ ,100.8
+ ,107.9
+ ,109.9
+ ,106.7
+ ,94.7
+ ,98.9
+ ,110.1
+ ,98.8
+ ,72.9
+ ,103.6
+ ,75.8
+ ,82.5
+ ,79.6
+ ,105.7
+ ,99.6
+ ,80.9
+ ,96.6
+ ,80
+ ,80.5
+ ,80.7
+ ,97.8
+ ,106.1
+ ,103.8
+ ,111.6
+ ,117.2
+ ,102.3
+ ,102.9
+ ,111.1
+ ,106.3
+ ,100.4
+ ,111.6
+ ,106.6
+ ,102.2
+ ,101.7
+ ,112
+ ,105.7
+ ,101.2
+ ,107
+ ,104.7
+ ,93.3
+ ,95.9
+ ,107.2
+ ,103.7
+ ,100.3
+ ,111.5
+ ,95.2
+ ,81
+ ,87.1
+ ,112.7
+ ,111.2
+ ,84.5
+ ,102
+ ,94
+ ,88.6
+ ,87.5
+ ,102.5
+ ,114.8
+ ,97.2
+ ,113.5
+ ,95.7
+ ,91.8
+ ,93.6
+ ,114.9
+ ,103.6
+ ,111.8
+ ,125.5
+ ,112.6
+ ,108.4
+ ,109.6
+ ,126.4
+ ,107
+ ,105.2
+ ,106.7
+ ,99.1
+ ,102.1
+ ,103.2
+ ,107.3
+ ,104.8
+ ,97.9
+ ,102.9
+ ,91.6
+ ,99.1
+ ,98.9
+ ,103.8
+ ,104.7
+ ,117.7
+ ,123.6
+ ,111.5
+ ,111.7
+ ,113.7
+ ,124.5
+ ,102
+ ,79.1
+ ,107.7
+ ,76.6
+ ,91.6
+ ,87.9
+ ,110.1
+ ,103.4
+ ,89.4
+ ,105.5
+ ,83.4
+ ,89.5
+ ,89.6
+ ,107.1
+ ,107
+ ,128.1
+ ,117.1
+ ,113.5
+ ,108.1
+ ,114.4
+ ,117.3
+ ,104
+ ,117.7
+ ,113.3
+ ,106.4
+ ,104.7
+ ,108.8
+ ,113.8
+ ,105.4
+ ,113.3
+ ,118
+ ,104.1
+ ,100.1
+ ,104.3
+ ,119
+ ,107.9
+ ,119.4
+ ,118.4
+ ,108.4
+ ,86.9
+ ,97
+ ,119.1
+ ,110.1
+ ,103.7
+ ,105.8
+ ,91
+ ,98.7
+ ,100.4
+ ,106.9
+ ,111
+ ,116
+ ,114.6
+ ,108.3
+ ,100.3
+ ,105.3
+ ,115
+ ,98.5
+ ,137.3
+ ,140.3
+ ,121
+ ,115.4
+ ,122.3
+ ,141.7
+ ,101.9
+ ,113.8
+ ,113.8
+ ,95.4
+ ,101.8
+ ,105.6
+ ,115.2
+ ,103.4
+ ,129.8
+ ,117.4
+ ,109.9
+ ,110.9
+ ,116.9
+ ,117.9
+ ,102.9
+ ,121.5
+ ,115.4
+ ,101.4
+ ,105.3
+ ,110.5
+ ,116.5
+ ,101
+ ,87
+ ,105.9
+ ,86
+ ,89.2
+ ,88.6
+ ,107.4
+ ,103.4
+ ,93.3
+ ,120.4
+ ,96.5
+ ,94.8
+ ,94.5
+ ,122.2
+ ,107.2
+ ,131.3
+ ,126.9
+ ,124.6
+ ,108.5
+ ,115.7
+ ,126.9
+ ,104.5
+ ,123.7
+ ,117.1
+ ,109.3
+ ,100.5
+ ,107.8
+ ,117.7
+ ,104.7
+ ,121.8
+ ,113.8
+ ,104.5
+ ,99.6
+ ,106.6
+ ,114.4
+ ,107
+ ,124.5
+ ,112.8
+ ,101.8
+ ,89.9
+ ,100.7
+ ,113.6
+ ,110.3
+ ,103.9
+ ,106.7
+ ,101.5
+ ,98.4
+ ,100.4
+ ,107
+ ,107.9
+ ,110.6
+ ,107.3
+ ,103.4
+ ,97.5
+ ,101.7
+ ,107.5
+ ,97.1
+ ,133.7
+ ,121.8
+ ,125.9
+ ,107
+ ,115.2
+ ,121.4
+ ,98.6
+ ,108
+ ,101.1
+ ,96.8
+ ,97.5
+ ,100.9
+ ,101.3
+ ,95.3
+ ,116.1
+ ,103.1
+ ,104.4
+ ,100.4
+ ,105.3
+ ,102.9
+ ,101.7
+ ,122.9
+ ,110.4
+ ,121.1
+ ,103.9
+ ,109.8
+ ,109.5
+ ,96.3
+ ,84.3
+ ,108.3
+ ,83.7
+ ,94.8
+ ,92.1
+ ,110.2
+ ,99
+ ,98.3
+ ,116.3
+ ,91.5
+ ,89.6
+ ,92.5
+ ,118.2
+ ,104)
+ ,dim=c(7
+ ,65)
+ ,dimnames=list(c('Investeringsgoederen'
+ ,'Consumptiegoederen'
+ ,'Duurzame-consumptiegoederen'
+ ,'Intermediaire-goederen'
+ ,'Intermediaire-en-investeringsgoederen'
+ ,'Niet-duurzame-consumptiegoederen'
+ ,'Energie')
+ ,1:65))
> y <- array(NA,dim=c(7,65),dimnames=list(c('Investeringsgoederen','Consumptiegoederen','Duurzame-consumptiegoederen','Intermediaire-goederen','Intermediaire-en-investeringsgoederen','Niet-duurzame-consumptiegoederen','Energie'),1:65))
> 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 = '3'
> 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
Duurzame-consumptiegoederen Investeringsgoederen Consumptiegoederen
1 95.0 105.3 97.2
2 101.4 110.0 100.9
3 112.8 129.2 109.4
4 76.4 99.9 101.3
5 93.5 100.6 101.2
6 117.6 117.2 104.8
7 118.6 136.0 111.4
8 110.1 126.0 106.9
9 113.1 125.0 103.9
10 91.2 112.5 107.6
11 106.2 119.3 104.7
12 115.5 121.1 109.1
13 106.2 127.1 109.3
14 95.9 116.0 102.2
15 113.2 131.9 109.8
16 78.3 101.0 106.2
17 79.8 92.6 95.1
18 121.2 127.2 118.7
19 125.6 124.3 116.9
20 97.2 103.8 105.3
21 102.8 106.4 119.5
22 88.8 84.2 96.5
23 95.3 91.9 99.3
24 107.6 103.4 113.8
25 95.0 93.0 102.7
26 87.5 110.6 98.8
27 106.7 107.9 109.9
28 75.8 72.9 103.6
29 80.0 80.9 96.6
30 117.2 103.8 111.6
31 106.6 100.4 111.6
32 104.7 101.2 107.0
33 95.2 100.3 111.5
34 94.0 84.5 102.0
35 95.7 97.2 113.5
36 112.6 111.8 125.5
37 99.1 105.2 106.7
38 91.6 97.9 102.9
39 111.5 117.7 123.6
40 76.6 79.1 107.7
41 83.4 89.4 105.5
42 113.5 128.1 117.1
43 106.4 117.7 113.3
44 104.1 113.3 118.0
45 108.4 119.4 118.4
46 91.0 103.7 105.8
47 108.3 116.0 114.6
48 121.0 137.3 140.3
49 95.4 113.8 113.8
50 109.9 129.8 117.4
51 101.4 121.5 115.4
52 86.0 87.0 105.9
53 96.5 93.3 120.4
54 124.6 131.3 126.9
55 109.3 123.7 117.1
56 104.5 121.8 113.8
57 101.8 124.5 112.8
58 101.5 103.9 106.7
59 103.4 110.6 107.3
60 125.9 133.7 121.8
61 96.8 108.0 101.1
62 104.4 116.1 103.1
63 121.1 122.9 110.4
64 83.7 84.3 108.3
65 91.5 98.3 116.3
Intermediaire-goederen Intermediaire-en-investeringsgoederen
1 109.3 108.0
2 115.7 113.8
3 118.5 122.1
4 100.9 100.6
5 99.9 100.2
6 109.9 112.3
7 120.7 125.7
8 113.3 117.5
9 96.4 105.8
10 116.2 115.1
11 119.2 119.3
12 117.9 119.0
13 126.5 126.8
14 116.3 116.3
15 124.7 127.1
16 108.8 106.5
17 96.7 95.5
18 118.4 121.3
19 115.4 118.3
20 91.7 95.6
21 82.5 90.1
22 81.7 82.6
23 84.5 86.9
24 92.4 95.9
25 86.8 88.8
26 85.1 93.2
27 94.7 98.9
28 82.5 79.6
29 80.5 80.7
30 102.3 102.9
31 102.2 101.7
32 93.3 95.9
33 81.0 87.1
34 88.6 87.5
35 91.8 93.6
36 108.4 109.6
37 102.1 103.2
38 99.1 98.9
39 111.7 113.7
40 91.6 87.9
41 89.5 89.6
42 108.1 114.4
43 104.7 108.8
44 100.1 104.3
45 86.9 97.0
46 98.7 100.4
47 100.3 105.3
48 115.4 122.3
49 101.8 105.6
50 110.9 116.9
51 105.3 110.5
52 89.2 88.6
53 94.8 94.5
54 108.5 115.7
55 100.5 107.8
56 99.6 106.6
57 89.9 100.7
58 98.4 100.4
59 97.5 101.7
60 107.0 115.2
61 97.5 100.9
62 100.4 105.3
63 103.9 109.8
64 94.8 92.1
65 89.6 92.5
Niet-duurzame-consumptiegoederen Energie
1 97.4 98.5
2 100.9 98.0
3 109.2 99.9
4 103.1 98.5
5 101.8 104.7
6 103.9 100.4
7 110.9 108.6
8 106.6 109.9
9 103.2 109.5
10 108.7 109.2
11 104.6 100.9
12 108.7 98.4
13 109.5 94.2
14 102.6 94.7
15 109.5 95.2
16 108.1 100.3
17 96.1 100.9
18 118.5 97.9
19 116.3 106.9
20 105.9 100.8
21 120.7 106.6
22 97.0 108.2
23 99.6 98.4
24 114.2 102.0
25 103.3 95.7
26 99.6 100.8
27 110.1 98.8
28 105.7 99.6
29 97.8 106.1
30 111.1 106.3
31 112.0 105.7
32 107.2 103.7
33 112.7 111.2
34 102.5 114.8
35 114.9 103.6
36 126.4 107.0
37 107.3 104.8
38 103.8 104.7
39 124.5 102.0
40 110.1 103.4
41 107.1 107.0
42 117.3 104.0
43 113.8 105.4
44 119.0 107.9
45 119.1 110.1
46 106.9 111.0
47 115.0 98.5
48 141.7 101.9
49 115.2 103.4
50 117.9 102.9
51 116.5 101.0
52 107.4 103.4
53 122.2 107.2
54 126.9 104.5
55 117.7 104.7
56 114.4 107.0
57 113.6 110.3
58 107.0 107.9
59 107.5 97.1
60 121.4 98.6
61 101.3 95.3
62 102.9 101.7
63 109.5 96.3
64 110.2 99.0
65 118.2 104.0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept)
0.12572
Investeringsgoederen
-0.48498
Consumptiegoederen
14.41963
`Intermediaire-goederen`
-1.01895
`Intermediaire-en-investeringsgoederen`
1.49938
`Niet-duurzame-consumptiegoederen`
-13.40901
Energie
-0.01146
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.36649 -0.49464 0.05398 0.66536 1.31847
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.12572 3.01172 0.042 0.967
Investeringsgoederen -0.48498 0.31948 -1.518 0.134
Consumptiegoederen 14.41963 0.21640 66.634 <2e-16
`Intermediaire-goederen` -1.01895 0.71015 -1.435 0.157
`Intermediaire-en-investeringsgoederen` 1.49938 1.02905 1.457 0.150
`Niet-duurzame-consumptiegoederen` -13.40901 0.20587 -65.132 <2e-16
Energie -0.01146 0.02484 -0.461 0.646
(Intercept)
Investeringsgoederen
Consumptiegoederen ***
`Intermediaire-goederen`
`Intermediaire-en-investeringsgoederen`
`Niet-duurzame-consumptiegoederen` ***
Energie
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8423 on 58 degrees of freedom
Multiple R-squared: 0.9961, Adjusted R-squared: 0.9957
F-statistic: 2484 on 6 and 58 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.4453991 0.89079823 0.55460088
[2,] 0.2876470 0.57529398 0.71235301
[3,] 0.2016569 0.40331379 0.79834311
[4,] 0.2228105 0.44562109 0.77718946
[5,] 0.1893085 0.37861694 0.81069153
[6,] 0.2350441 0.47008812 0.76495594
[7,] 0.3335832 0.66716648 0.66641676
[8,] 0.4525356 0.90507111 0.54746444
[9,] 0.5853833 0.82923332 0.41461666
[10,] 0.6767380 0.64652399 0.32326200
[11,] 0.6067746 0.78645070 0.39322535
[12,] 0.5615433 0.87691345 0.43845673
[13,] 0.5122645 0.97547110 0.48773555
[14,] 0.4467600 0.89351991 0.55324005
[15,] 0.4738843 0.94776857 0.52611572
[16,] 0.4719473 0.94389453 0.52805273
[17,] 0.4379842 0.87596842 0.56201579
[18,] 0.4336641 0.86732820 0.56633590
[19,] 0.7573766 0.48524684 0.24262342
[20,] 0.7287177 0.54256457 0.27128228
[21,] 0.8082276 0.38354489 0.19177245
[22,] 0.8406822 0.31863556 0.15931778
[23,] 0.9130280 0.17394401 0.08697201
[24,] 0.8849538 0.23009248 0.11504624
[25,] 0.8568093 0.28638135 0.14319068
[26,] 0.9391792 0.12164161 0.06082080
[27,] 0.9163522 0.16729562 0.08364781
[28,] 0.9381261 0.12374770 0.06187385
[29,] 0.9741853 0.05162934 0.02581467
[30,] 0.9654360 0.06912795 0.03456398
[31,] 0.9738515 0.05229706 0.02614853
[32,] 0.9822183 0.03556344 0.01778172
[33,] 0.9715912 0.05681750 0.02840875
[34,] 0.9651795 0.06964093 0.03482046
[35,] 0.9425279 0.11494412 0.05747206
[36,] 0.9086912 0.18261763 0.09130882
[37,] 0.8707631 0.25847379 0.12923690
[38,] 0.8269526 0.34609479 0.17304740
[39,] 0.8041391 0.39172173 0.19586087
[40,] 0.7412227 0.51755455 0.25877727
[41,] 0.6595888 0.68082240 0.34041120
[42,] 0.6681464 0.66370727 0.33185363
[43,] 0.5506878 0.89862432 0.44931216
[44,] 0.4255532 0.85110631 0.57444685
[45,] 0.7212294 0.55754120 0.27877060
[46,] 0.7690659 0.46186824 0.23093412
> postscript(file="/var/wessaorg/rcomp/tmp/17pzz1353442514.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/wessaorg/rcomp/tmp/2worp1353442514.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/wessaorg/rcomp/tmp/31tz01353442514.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/wessaorg/rcomp/tmp/4tkno1353442514.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/wessaorg/rcomp/tmp/5gmbm1353442514.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 = 65
Frequency = 1
1 2 3 4 5 6
0.95895894 1.03640677 0.90587761 -0.41282113 0.68876865 1.08537899
7 8 9 10 11 12
0.90340304 -0.44719828 0.05397937 -1.28376746 0.51840726 1.31846654
13 14 15 16 17 18
-0.20871152 -0.67883128 -0.36309956 -2.36648885 -1.61976780 0.01184995
19 20 21 22 23 24
1.00535499 0.29421582 -0.21152872 -0.67154619 0.34466411 -0.49464082
25 26 27 28 29 30
0.62855681 0.01632860 -0.14385282 0.34196920 -0.18471211 -0.90416608
31 32 33 34 35 36
0.60550260 0.66535901 0.33749767 -1.12503189 1.16623185 -0.72158128
37 38 39 40 41 42
0.70595168 0.91801183 0.11774165 0.89912981 -0.45693403 -0.35064017
43 44 45 46 47 48
0.31679443 -0.07382792 0.27812246 0.29742897 -0.57699693 -0.17525004
49 50 51 52 53 54
0.80831943 -0.31465766 1.09481103 0.38641267 0.21381119 -1.82047208
55 56 57 58 59 60
0.83896900 -0.63888470 0.66331410 -1.08354895 -0.87148748 -0.41231536
61 62 63 64 65
-1.28770910 -0.71320945 -0.72195858 0.12295729 0.78668494
> postscript(file="/var/wessaorg/rcomp/tmp/6g4p91353442514.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 = 65
Frequency = 1
lag(myerror, k = 1) myerror
0 0.95895894 NA
1 1.03640677 0.95895894
2 0.90587761 1.03640677
3 -0.41282113 0.90587761
4 0.68876865 -0.41282113
5 1.08537899 0.68876865
6 0.90340304 1.08537899
7 -0.44719828 0.90340304
8 0.05397937 -0.44719828
9 -1.28376746 0.05397937
10 0.51840726 -1.28376746
11 1.31846654 0.51840726
12 -0.20871152 1.31846654
13 -0.67883128 -0.20871152
14 -0.36309956 -0.67883128
15 -2.36648885 -0.36309956
16 -1.61976780 -2.36648885
17 0.01184995 -1.61976780
18 1.00535499 0.01184995
19 0.29421582 1.00535499
20 -0.21152872 0.29421582
21 -0.67154619 -0.21152872
22 0.34466411 -0.67154619
23 -0.49464082 0.34466411
24 0.62855681 -0.49464082
25 0.01632860 0.62855681
26 -0.14385282 0.01632860
27 0.34196920 -0.14385282
28 -0.18471211 0.34196920
29 -0.90416608 -0.18471211
30 0.60550260 -0.90416608
31 0.66535901 0.60550260
32 0.33749767 0.66535901
33 -1.12503189 0.33749767
34 1.16623185 -1.12503189
35 -0.72158128 1.16623185
36 0.70595168 -0.72158128
37 0.91801183 0.70595168
38 0.11774165 0.91801183
39 0.89912981 0.11774165
40 -0.45693403 0.89912981
41 -0.35064017 -0.45693403
42 0.31679443 -0.35064017
43 -0.07382792 0.31679443
44 0.27812246 -0.07382792
45 0.29742897 0.27812246
46 -0.57699693 0.29742897
47 -0.17525004 -0.57699693
48 0.80831943 -0.17525004
49 -0.31465766 0.80831943
50 1.09481103 -0.31465766
51 0.38641267 1.09481103
52 0.21381119 0.38641267
53 -1.82047208 0.21381119
54 0.83896900 -1.82047208
55 -0.63888470 0.83896900
56 0.66331410 -0.63888470
57 -1.08354895 0.66331410
58 -0.87148748 -1.08354895
59 -0.41231536 -0.87148748
60 -1.28770910 -0.41231536
61 -0.71320945 -1.28770910
62 -0.72195858 -0.71320945
63 0.12295729 -0.72195858
64 0.78668494 0.12295729
65 NA 0.78668494
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1.03640677 0.95895894
[2,] 0.90587761 1.03640677
[3,] -0.41282113 0.90587761
[4,] 0.68876865 -0.41282113
[5,] 1.08537899 0.68876865
[6,] 0.90340304 1.08537899
[7,] -0.44719828 0.90340304
[8,] 0.05397937 -0.44719828
[9,] -1.28376746 0.05397937
[10,] 0.51840726 -1.28376746
[11,] 1.31846654 0.51840726
[12,] -0.20871152 1.31846654
[13,] -0.67883128 -0.20871152
[14,] -0.36309956 -0.67883128
[15,] -2.36648885 -0.36309956
[16,] -1.61976780 -2.36648885
[17,] 0.01184995 -1.61976780
[18,] 1.00535499 0.01184995
[19,] 0.29421582 1.00535499
[20,] -0.21152872 0.29421582
[21,] -0.67154619 -0.21152872
[22,] 0.34466411 -0.67154619
[23,] -0.49464082 0.34466411
[24,] 0.62855681 -0.49464082
[25,] 0.01632860 0.62855681
[26,] -0.14385282 0.01632860
[27,] 0.34196920 -0.14385282
[28,] -0.18471211 0.34196920
[29,] -0.90416608 -0.18471211
[30,] 0.60550260 -0.90416608
[31,] 0.66535901 0.60550260
[32,] 0.33749767 0.66535901
[33,] -1.12503189 0.33749767
[34,] 1.16623185 -1.12503189
[35,] -0.72158128 1.16623185
[36,] 0.70595168 -0.72158128
[37,] 0.91801183 0.70595168
[38,] 0.11774165 0.91801183
[39,] 0.89912981 0.11774165
[40,] -0.45693403 0.89912981
[41,] -0.35064017 -0.45693403
[42,] 0.31679443 -0.35064017
[43,] -0.07382792 0.31679443
[44,] 0.27812246 -0.07382792
[45,] 0.29742897 0.27812246
[46,] -0.57699693 0.29742897
[47,] -0.17525004 -0.57699693
[48,] 0.80831943 -0.17525004
[49,] -0.31465766 0.80831943
[50,] 1.09481103 -0.31465766
[51,] 0.38641267 1.09481103
[52,] 0.21381119 0.38641267
[53,] -1.82047208 0.21381119
[54,] 0.83896900 -1.82047208
[55,] -0.63888470 0.83896900
[56,] 0.66331410 -0.63888470
[57,] -1.08354895 0.66331410
[58,] -0.87148748 -1.08354895
[59,] -0.41231536 -0.87148748
[60,] -1.28770910 -0.41231536
[61,] -0.71320945 -1.28770910
[62,] -0.72195858 -0.71320945
[63,] 0.12295729 -0.72195858
[64,] 0.78668494 0.12295729
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1.03640677 0.95895894
2 0.90587761 1.03640677
3 -0.41282113 0.90587761
4 0.68876865 -0.41282113
5 1.08537899 0.68876865
6 0.90340304 1.08537899
7 -0.44719828 0.90340304
8 0.05397937 -0.44719828
9 -1.28376746 0.05397937
10 0.51840726 -1.28376746
11 1.31846654 0.51840726
12 -0.20871152 1.31846654
13 -0.67883128 -0.20871152
14 -0.36309956 -0.67883128
15 -2.36648885 -0.36309956
16 -1.61976780 -2.36648885
17 0.01184995 -1.61976780
18 1.00535499 0.01184995
19 0.29421582 1.00535499
20 -0.21152872 0.29421582
21 -0.67154619 -0.21152872
22 0.34466411 -0.67154619
23 -0.49464082 0.34466411
24 0.62855681 -0.49464082
25 0.01632860 0.62855681
26 -0.14385282 0.01632860
27 0.34196920 -0.14385282
28 -0.18471211 0.34196920
29 -0.90416608 -0.18471211
30 0.60550260 -0.90416608
31 0.66535901 0.60550260
32 0.33749767 0.66535901
33 -1.12503189 0.33749767
34 1.16623185 -1.12503189
35 -0.72158128 1.16623185
36 0.70595168 -0.72158128
37 0.91801183 0.70595168
38 0.11774165 0.91801183
39 0.89912981 0.11774165
40 -0.45693403 0.89912981
41 -0.35064017 -0.45693403
42 0.31679443 -0.35064017
43 -0.07382792 0.31679443
44 0.27812246 -0.07382792
45 0.29742897 0.27812246
46 -0.57699693 0.29742897
47 -0.17525004 -0.57699693
48 0.80831943 -0.17525004
49 -0.31465766 0.80831943
50 1.09481103 -0.31465766
51 0.38641267 1.09481103
52 0.21381119 0.38641267
53 -1.82047208 0.21381119
54 0.83896900 -1.82047208
55 -0.63888470 0.83896900
56 0.66331410 -0.63888470
57 -1.08354895 0.66331410
58 -0.87148748 -1.08354895
59 -0.41231536 -0.87148748
60 -1.28770910 -0.41231536
61 -0.71320945 -1.28770910
62 -0.72195858 -0.71320945
63 0.12295729 -0.72195858
64 0.78668494 0.12295729
> 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/wessaorg/rcomp/tmp/7ge041353442514.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/wessaorg/rcomp/tmp/8emuw1353442514.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/wessaorg/rcomp/tmp/9d8nx1353442514.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/wessaorg/rcomp/tmp/10d04s1353442514.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11j42o1353442514.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/wessaorg/rcomp/tmp/12w40v1353442514.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/wessaorg/rcomp/tmp/138kh71353442514.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/wessaorg/rcomp/tmp/14izbm1353442514.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/wessaorg/rcomp/tmp/15ege01353442514.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/wessaorg/rcomp/tmp/1661hq1353442514.tab")
+ }
>
> try(system("convert tmp/17pzz1353442514.ps tmp/17pzz1353442514.png",intern=TRUE))
character(0)
> try(system("convert tmp/2worp1353442514.ps tmp/2worp1353442514.png",intern=TRUE))
character(0)
> try(system("convert tmp/31tz01353442514.ps tmp/31tz01353442514.png",intern=TRUE))
character(0)
> try(system("convert tmp/4tkno1353442514.ps tmp/4tkno1353442514.png",intern=TRUE))
character(0)
> try(system("convert tmp/5gmbm1353442514.ps tmp/5gmbm1353442514.png",intern=TRUE))
character(0)
> try(system("convert tmp/6g4p91353442514.ps tmp/6g4p91353442514.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ge041353442514.ps tmp/7ge041353442514.png",intern=TRUE))
character(0)
> try(system("convert tmp/8emuw1353442514.ps tmp/8emuw1353442514.png",intern=TRUE))
character(0)
> try(system("convert tmp/9d8nx1353442514.ps tmp/9d8nx1353442514.png",intern=TRUE))
character(0)
> try(system("convert tmp/10d04s1353442514.ps tmp/10d04s1353442514.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
7.594 1.417 8.966