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(101.9
+ ,122.2
+ ,19
+ ,73
+ ,77.8
+ ,74.8
+ ,80.2
+ ,102
+ ,123.7
+ ,22
+ ,72
+ ,73
+ ,77.8
+ ,74.8
+ ,100.7
+ ,122.6
+ ,23
+ ,75.8
+ ,72
+ ,73
+ ,77.8
+ ,99
+ ,115.7
+ ,20
+ ,72.6
+ ,75.8
+ ,72
+ ,73
+ ,96.5
+ ,116.1
+ ,14
+ ,71.9
+ ,72.6
+ ,75.8
+ ,72
+ ,101.8
+ ,120.5
+ ,14
+ ,74.8
+ ,71.9
+ ,72.6
+ ,75.8
+ ,100.5
+ ,122.6
+ ,14
+ ,72.9
+ ,74.8
+ ,71.9
+ ,72.6
+ ,103.3
+ ,119.9
+ ,15
+ ,72.9
+ ,72.9
+ ,74.8
+ ,71.9
+ ,102.3
+ ,120.7
+ ,11
+ ,79.9
+ ,72.9
+ ,72.9
+ ,74.8
+ ,100.4
+ ,120.2
+ ,17
+ ,74
+ ,79.9
+ ,72.9
+ ,72.9
+ ,103
+ ,122.1
+ ,16
+ ,76
+ ,74
+ ,79.9
+ ,72.9
+ ,99
+ ,119.3
+ ,20
+ ,69.6
+ ,76
+ ,74
+ ,79.9
+ ,104.8
+ ,121.7
+ ,24
+ ,77.3
+ ,69.6
+ ,76
+ ,74
+ ,104.5
+ ,113.5
+ ,23
+ ,75.2
+ ,77.3
+ ,69.6
+ ,76
+ ,104.8
+ ,123.7
+ ,20
+ ,75.8
+ ,75.2
+ ,77.3
+ ,69.6
+ ,103.8
+ ,123.4
+ ,21
+ ,77.6
+ ,75.8
+ ,75.2
+ ,77.3
+ ,106.3
+ ,126.4
+ ,19
+ ,76.7
+ ,77.6
+ ,75.8
+ ,75.2
+ ,105.2
+ ,124.1
+ ,23
+ ,77
+ ,76.7
+ ,77.6
+ ,75.8
+ ,108.2
+ ,125.6
+ ,23
+ ,77.9
+ ,77
+ ,76.7
+ ,77.6
+ ,106.2
+ ,124.8
+ ,23
+ ,76.7
+ ,77.9
+ ,77
+ ,76.7
+ ,103.9
+ ,123
+ ,23
+ ,71.9
+ ,76.7
+ ,77.9
+ ,77
+ ,104.9
+ ,126.9
+ ,27
+ ,73.4
+ ,71.9
+ ,76.7
+ ,77.9
+ ,106.2
+ ,127.3
+ ,26
+ ,72.5
+ ,73.4
+ ,71.9
+ ,76.7
+ ,107.9
+ ,129
+ ,17
+ ,73.7
+ ,72.5
+ ,73.4
+ ,71.9
+ ,106.9
+ ,126.2
+ ,24
+ ,69.5
+ ,73.7
+ ,72.5
+ ,73.4
+ ,110.3
+ ,125.4
+ ,26
+ ,74.7
+ ,69.5
+ ,73.7
+ ,72.5
+ ,109.8
+ ,126.3
+ ,24
+ ,72.5
+ ,74.7
+ ,69.5
+ ,73.7
+ ,108.3
+ ,126.3
+ ,27
+ ,72.1
+ ,72.5
+ ,74.7
+ ,69.5
+ ,110.9
+ ,128.4
+ ,27
+ ,70.7
+ ,72.1
+ ,72.5
+ ,74.7
+ ,109.8
+ ,127.2
+ ,26
+ ,71.4
+ ,70.7
+ ,72.1
+ ,72.5
+ ,109.3
+ ,128.5
+ ,24
+ ,69.5
+ ,71.4
+ ,70.7
+ ,72.1
+ ,109
+ ,129
+ ,23
+ ,73.5
+ ,69.5
+ ,71.4
+ ,70.7
+ ,107.9
+ ,128.9
+ ,23
+ ,72.4
+ ,73.5
+ ,69.5
+ ,71.4
+ ,108.4
+ ,128.3
+ ,24
+ ,74.5
+ ,72.4
+ ,73.5
+ ,69.5
+ ,107.2
+ ,124.6
+ ,17
+ ,72.2
+ ,74.5
+ ,72.4
+ ,73.5
+ ,109.5
+ ,126.2
+ ,21
+ ,73
+ ,72.2
+ ,74.5
+ ,72.4
+ ,109.9
+ ,129.1
+ ,19
+ ,73.3
+ ,73
+ ,72.2
+ ,74.5
+ ,108
+ ,127.3
+ ,22
+ ,71.3
+ ,73.3
+ ,73
+ ,72.2
+ ,114.7
+ ,129.2
+ ,22
+ ,73.6
+ ,71.3
+ ,73.3
+ ,73
+ ,115.6
+ ,130.4
+ ,18
+ ,71.3
+ ,73.6
+ ,71.3
+ ,73.3
+ ,107.6
+ ,125.9
+ ,16
+ ,71.2
+ ,71.3
+ ,73.6
+ ,71.3
+ ,115.9
+ ,135.8
+ ,14
+ ,81.4
+ ,71.2
+ ,71.3
+ ,73.6
+ ,111.8
+ ,126.4
+ ,12
+ ,76.1
+ ,81.4
+ ,71.2
+ ,71.3
+ ,110
+ ,129.5
+ ,14
+ ,71.1
+ ,76.1
+ ,81.4
+ ,71.2
+ ,109.2
+ ,128.4
+ ,16
+ ,75.7
+ ,71.1
+ ,76.1
+ ,81.4
+ ,108
+ ,125.6
+ ,8
+ ,70
+ ,75.7
+ ,71.1
+ ,76.1
+ ,105.6
+ ,127.7
+ ,3
+ ,68.5
+ ,70
+ ,75.7
+ ,71.1
+ ,103
+ ,126.4
+ ,0
+ ,56.7
+ ,68.5
+ ,70
+ ,75.7
+ ,99.6
+ ,124.2
+ ,5
+ ,57.9
+ ,56.7
+ ,68.5
+ ,70
+ ,97.9
+ ,126.4
+ ,1
+ ,58.8
+ ,57.9
+ ,56.7
+ ,68.5
+ ,97.6
+ ,123.7
+ ,1
+ ,59.3
+ ,58.8
+ ,57.9
+ ,56.7
+ ,96.2
+ ,121.8
+ ,3
+ ,61.3
+ ,59.3
+ ,58.8
+ ,57.9
+ ,97.9
+ ,124
+ ,6
+ ,62.9
+ ,61.3
+ ,59.3
+ ,58.8
+ ,94.5
+ ,122.7
+ ,7
+ ,61.4
+ ,62.9
+ ,61.3
+ ,59.3
+ ,95.4
+ ,122.9
+ ,8
+ ,64.5
+ ,61.4
+ ,62.9
+ ,61.3
+ ,94.4
+ ,121
+ ,14
+ ,63.8
+ ,64.5
+ ,61.4
+ ,62.9
+ ,96.3
+ ,122.8
+ ,14
+ ,61.6
+ ,63.8
+ ,64.5
+ ,61.4
+ ,95.1
+ ,122.9
+ ,13
+ ,64.7
+ ,61.6
+ ,63.8
+ ,64.5)
+ ,dim=c(7
+ ,58)
+ ,dimnames=list(c('totid'
+ ,'ndzcg'
+ ,'indc'
+ ,'Y'
+ ,'y1'
+ ,'y2'
+ ,'y3
')
+ ,1:58))
> y <- array(NA,dim=c(7,58),dimnames=list(c('totid','ndzcg','indc','Y','y1','y2','y3
'),1:58))
> 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 = '4'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Y totid ndzcg indc y1 y2 y3\r M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 73.0 101.9 122.2 19 77.8 74.8 80.2 1 0 0 0 0 0 0 0 0 0 0 1
2 72.0 102.0 123.7 22 73.0 77.8 74.8 0 1 0 0 0 0 0 0 0 0 0 2
3 75.8 100.7 122.6 23 72.0 73.0 77.8 0 0 1 0 0 0 0 0 0 0 0 3
4 72.6 99.0 115.7 20 75.8 72.0 73.0 0 0 0 1 0 0 0 0 0 0 0 4
5 71.9 96.5 116.1 14 72.6 75.8 72.0 0 0 0 0 1 0 0 0 0 0 0 5
6 74.8 101.8 120.5 14 71.9 72.6 75.8 0 0 0 0 0 1 0 0 0 0 0 6
7 72.9 100.5 122.6 14 74.8 71.9 72.6 0 0 0 0 0 0 1 0 0 0 0 7
8 72.9 103.3 119.9 15 72.9 74.8 71.9 0 0 0 0 0 0 0 1 0 0 0 8
9 79.9 102.3 120.7 11 72.9 72.9 74.8 0 0 0 0 0 0 0 0 1 0 0 9
10 74.0 100.4 120.2 17 79.9 72.9 72.9 0 0 0 0 0 0 0 0 0 1 0 10
11 76.0 103.0 122.1 16 74.0 79.9 72.9 0 0 0 0 0 0 0 0 0 0 1 11
12 69.6 99.0 119.3 20 76.0 74.0 79.9 0 0 0 0 0 0 0 0 0 0 0 12
13 77.3 104.8 121.7 24 69.6 76.0 74.0 1 0 0 0 0 0 0 0 0 0 0 13
14 75.2 104.5 113.5 23 77.3 69.6 76.0 0 1 0 0 0 0 0 0 0 0 0 14
15 75.8 104.8 123.7 20 75.2 77.3 69.6 0 0 1 0 0 0 0 0 0 0 0 15
16 77.6 103.8 123.4 21 75.8 75.2 77.3 0 0 0 1 0 0 0 0 0 0 0 16
17 76.7 106.3 126.4 19 77.6 75.8 75.2 0 0 0 0 1 0 0 0 0 0 0 17
18 77.0 105.2 124.1 23 76.7 77.6 75.8 0 0 0 0 0 1 0 0 0 0 0 18
19 77.9 108.2 125.6 23 77.0 76.7 77.6 0 0 0 0 0 0 1 0 0 0 0 19
20 76.7 106.2 124.8 23 77.9 77.0 76.7 0 0 0 0 0 0 0 1 0 0 0 20
21 71.9 103.9 123.0 23 76.7 77.9 77.0 0 0 0 0 0 0 0 0 1 0 0 21
22 73.4 104.9 126.9 27 71.9 76.7 77.9 0 0 0 0 0 0 0 0 0 1 0 22
23 72.5 106.2 127.3 26 73.4 71.9 76.7 0 0 0 0 0 0 0 0 0 0 1 23
24 73.7 107.9 129.0 17 72.5 73.4 71.9 0 0 0 0 0 0 0 0 0 0 0 24
25 69.5 106.9 126.2 24 73.7 72.5 73.4 1 0 0 0 0 0 0 0 0 0 0 25
26 74.7 110.3 125.4 26 69.5 73.7 72.5 0 1 0 0 0 0 0 0 0 0 0 26
27 72.5 109.8 126.3 24 74.7 69.5 73.7 0 0 1 0 0 0 0 0 0 0 0 27
28 72.1 108.3 126.3 27 72.5 74.7 69.5 0 0 0 1 0 0 0 0 0 0 0 28
29 70.7 110.9 128.4 27 72.1 72.5 74.7 0 0 0 0 1 0 0 0 0 0 0 29
30 71.4 109.8 127.2 26 70.7 72.1 72.5 0 0 0 0 0 1 0 0 0 0 0 30
31 69.5 109.3 128.5 24 71.4 70.7 72.1 0 0 0 0 0 0 1 0 0 0 0 31
32 73.5 109.0 129.0 23 69.5 71.4 70.7 0 0 0 0 0 0 0 1 0 0 0 32
33 72.4 107.9 128.9 23 73.5 69.5 71.4 0 0 0 0 0 0 0 0 1 0 0 33
34 74.5 108.4 128.3 24 72.4 73.5 69.5 0 0 0 0 0 0 0 0 0 1 0 34
35 72.2 107.2 124.6 17 74.5 72.4 73.5 0 0 0 0 0 0 0 0 0 0 1 35
36 73.0 109.5 126.2 21 72.2 74.5 72.4 0 0 0 0 0 0 0 0 0 0 0 36
37 73.3 109.9 129.1 19 73.0 72.2 74.5 1 0 0 0 0 0 0 0 0 0 0 37
38 71.3 108.0 127.3 22 73.3 73.0 72.2 0 1 0 0 0 0 0 0 0 0 0 38
39 73.6 114.7 129.2 22 71.3 73.3 73.0 0 0 1 0 0 0 0 0 0 0 0 39
40 71.3 115.6 130.4 18 73.6 71.3 73.3 0 0 0 1 0 0 0 0 0 0 0 40
41 71.2 107.6 125.9 16 71.3 73.6 71.3 0 0 0 0 1 0 0 0 0 0 0 41
42 81.4 115.9 135.8 14 71.2 71.3 73.6 0 0 0 0 0 1 0 0 0 0 0 42
43 76.1 111.8 126.4 12 81.4 71.2 71.3 0 0 0 0 0 0 1 0 0 0 0 43
44 71.1 110.0 129.5 14 76.1 81.4 71.2 0 0 0 0 0 0 0 1 0 0 0 44
45 75.7 109.2 128.4 16 71.1 76.1 81.4 0 0 0 0 0 0 0 0 1 0 0 45
46 70.0 108.0 125.6 8 75.7 71.1 76.1 0 0 0 0 0 0 0 0 0 1 0 46
47 68.5 105.6 127.7 3 70.0 75.7 71.1 0 0 0 0 0 0 0 0 0 0 1 47
48 56.7 103.0 126.4 0 68.5 70.0 75.7 0 0 0 0 0 0 0 0 0 0 0 48
49 57.9 99.6 124.2 5 56.7 68.5 70.0 1 0 0 0 0 0 0 0 0 0 0 49
50 58.8 97.9 126.4 1 57.9 56.7 68.5 0 1 0 0 0 0 0 0 0 0 0 50
51 59.3 97.6 123.7 1 58.8 57.9 56.7 0 0 1 0 0 0 0 0 0 0 0 51
52 61.3 96.2 121.8 3 59.3 58.8 57.9 0 0 0 1 0 0 0 0 0 0 0 52
53 62.9 97.9 124.0 6 61.3 59.3 58.8 0 0 0 0 1 0 0 0 0 0 0 53
54 61.4 94.5 122.7 7 62.9 61.3 59.3 0 0 0 0 0 1 0 0 0 0 0 54
55 64.5 95.4 122.9 8 61.4 62.9 61.3 0 0 0 0 0 0 1 0 0 0 0 55
56 63.8 94.4 121.0 14 64.5 61.4 62.9 0 0 0 0 0 0 0 1 0 0 0 56
57 61.6 96.3 122.8 14 63.8 64.5 61.4 0 0 0 0 0 0 0 0 1 0 0 57
58 64.7 95.1 122.9 13 61.6 63.8 64.5 0 0 0 0 0 0 0 0 0 1 0 58
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) totid ndzcg indc y1 y2
19.54471 0.47415 -0.08721 0.04613 0.16357 0.11584
`y3\r` M1 M2 M3 M4 M5
-0.08371 1.37117 1.77995 2.41302 2.28548 2.58506
M6 M7 M8 M9 M10 M11
4.83378 3.70224 3.24367 4.85075 4.14558 3.09305
t
-0.15870
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-6.01795 -1.52868 -0.03425 1.32187 6.18539
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.54471 18.44658 1.060 0.29588
totid 0.47415 0.18445 2.571 0.01408 *
ndzcg -0.08721 0.21140 -0.413 0.68222
indc 0.04613 0.07952 0.580 0.56519
y1 0.16357 0.15822 1.034 0.30762
y2 0.11584 0.14870 0.779 0.44070
`y3\r` -0.08371 0.14405 -0.581 0.56451
M1 1.37117 2.05870 0.666 0.50931
M2 1.77995 2.13776 0.833 0.41013
M3 2.41302 2.17547 1.109 0.27414
M4 2.28548 2.14510 1.065 0.29323
M5 2.58506 2.08146 1.242 0.22167
M6 4.83378 2.07254 2.332 0.02494 *
M7 3.70224 2.09409 1.768 0.08489 .
M8 3.24367 2.09148 1.551 0.12900
M9 4.85075 2.02260 2.398 0.02135 *
M10 4.14558 2.06055 2.012 0.05118 .
M11 3.09305 2.13337 1.450 0.15510
t -0.15870 0.05229 -3.035 0.00427 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.965 on 39 degrees of freedom
Multiple R-squared: 0.8108, Adjusted R-squared: 0.7234
F-statistic: 9.282 on 18 and 39 DF, p-value: 4.307e-09
> 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.6658540 0.6682920 0.3341460
[2,] 0.5342101 0.9315798 0.4657899
[3,] 0.4631266 0.9262532 0.5368734
[4,] 0.4220059 0.8440119 0.5779941
[5,] 0.3149054 0.6298109 0.6850946
[6,] 0.3348237 0.6696473 0.6651763
[7,] 0.2834426 0.5668852 0.7165574
[8,] 0.2513870 0.5027739 0.7486130
[9,] 0.1957585 0.3915171 0.8042415
[10,] 0.2483889 0.4967778 0.7516111
[11,] 0.2012342 0.4024684 0.7987658
[12,] 0.1594257 0.3188514 0.8405743
[13,] 0.1747520 0.3495039 0.8252480
[14,] 0.3857181 0.7714362 0.6142819
[15,] 0.2740704 0.5481408 0.7259296
> postscript(file="/var/www/html/rcomp/tmp/1vn9v1258660722.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/26r0t1258660722.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/3xcq11258660722.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/41n051258660722.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/5geq91258660722.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 = 58
Frequency = 1
1 2 3 4 5 6
-0.96985742 -2.28934643 2.48133742 -0.99720443 -0.34150497 -0.85757508
7 8 9 10 11 12
-1.32890930 -2.40460004 4.33829782 -1.42131796 0.92312783 0.18504947
13 14 15 16 17 18
4.26853933 1.04105413 0.96816803 4.24593307 1.83385285 0.16923357
19 20 21 22 23 24
1.27368677 1.31220692 -3.88545304 -0.84069324 -0.85462292 2.92614376
25 26 27 28 29 30
-2.54572057 1.10268555 -1.42743034 -1.56242336 -3.39741667 -4.23321850
31 32 33 34 35 36
-4.38607261 0.57567878 -1.83542970 0.35044691 -0.05032172 2.90675230
37 38 39 40 41 42
2.46113885 0.48232510 -0.34382746 -2.61457472 0.57997788 6.18538957
43 44 45 46 47 48
1.54287801 -2.13129313 3.49703496 -1.26211579 -0.01818319 -6.01794552
49 50 51 52 53 54
-3.21410020 -0.33671835 -1.67824766 0.92826945 1.32509091 -1.26382955
55 56 57 58
2.89841712 2.64800747 -2.11445003 3.17368008
> postscript(file="/var/www/html/rcomp/tmp/67wuo1258660722.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 = 58
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.96985742 NA
1 -2.28934643 -0.96985742
2 2.48133742 -2.28934643
3 -0.99720443 2.48133742
4 -0.34150497 -0.99720443
5 -0.85757508 -0.34150497
6 -1.32890930 -0.85757508
7 -2.40460004 -1.32890930
8 4.33829782 -2.40460004
9 -1.42131796 4.33829782
10 0.92312783 -1.42131796
11 0.18504947 0.92312783
12 4.26853933 0.18504947
13 1.04105413 4.26853933
14 0.96816803 1.04105413
15 4.24593307 0.96816803
16 1.83385285 4.24593307
17 0.16923357 1.83385285
18 1.27368677 0.16923357
19 1.31220692 1.27368677
20 -3.88545304 1.31220692
21 -0.84069324 -3.88545304
22 -0.85462292 -0.84069324
23 2.92614376 -0.85462292
24 -2.54572057 2.92614376
25 1.10268555 -2.54572057
26 -1.42743034 1.10268555
27 -1.56242336 -1.42743034
28 -3.39741667 -1.56242336
29 -4.23321850 -3.39741667
30 -4.38607261 -4.23321850
31 0.57567878 -4.38607261
32 -1.83542970 0.57567878
33 0.35044691 -1.83542970
34 -0.05032172 0.35044691
35 2.90675230 -0.05032172
36 2.46113885 2.90675230
37 0.48232510 2.46113885
38 -0.34382746 0.48232510
39 -2.61457472 -0.34382746
40 0.57997788 -2.61457472
41 6.18538957 0.57997788
42 1.54287801 6.18538957
43 -2.13129313 1.54287801
44 3.49703496 -2.13129313
45 -1.26211579 3.49703496
46 -0.01818319 -1.26211579
47 -6.01794552 -0.01818319
48 -3.21410020 -6.01794552
49 -0.33671835 -3.21410020
50 -1.67824766 -0.33671835
51 0.92826945 -1.67824766
52 1.32509091 0.92826945
53 -1.26382955 1.32509091
54 2.89841712 -1.26382955
55 2.64800747 2.89841712
56 -2.11445003 2.64800747
57 3.17368008 -2.11445003
58 NA 3.17368008
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2.28934643 -0.96985742
[2,] 2.48133742 -2.28934643
[3,] -0.99720443 2.48133742
[4,] -0.34150497 -0.99720443
[5,] -0.85757508 -0.34150497
[6,] -1.32890930 -0.85757508
[7,] -2.40460004 -1.32890930
[8,] 4.33829782 -2.40460004
[9,] -1.42131796 4.33829782
[10,] 0.92312783 -1.42131796
[11,] 0.18504947 0.92312783
[12,] 4.26853933 0.18504947
[13,] 1.04105413 4.26853933
[14,] 0.96816803 1.04105413
[15,] 4.24593307 0.96816803
[16,] 1.83385285 4.24593307
[17,] 0.16923357 1.83385285
[18,] 1.27368677 0.16923357
[19,] 1.31220692 1.27368677
[20,] -3.88545304 1.31220692
[21,] -0.84069324 -3.88545304
[22,] -0.85462292 -0.84069324
[23,] 2.92614376 -0.85462292
[24,] -2.54572057 2.92614376
[25,] 1.10268555 -2.54572057
[26,] -1.42743034 1.10268555
[27,] -1.56242336 -1.42743034
[28,] -3.39741667 -1.56242336
[29,] -4.23321850 -3.39741667
[30,] -4.38607261 -4.23321850
[31,] 0.57567878 -4.38607261
[32,] -1.83542970 0.57567878
[33,] 0.35044691 -1.83542970
[34,] -0.05032172 0.35044691
[35,] 2.90675230 -0.05032172
[36,] 2.46113885 2.90675230
[37,] 0.48232510 2.46113885
[38,] -0.34382746 0.48232510
[39,] -2.61457472 -0.34382746
[40,] 0.57997788 -2.61457472
[41,] 6.18538957 0.57997788
[42,] 1.54287801 6.18538957
[43,] -2.13129313 1.54287801
[44,] 3.49703496 -2.13129313
[45,] -1.26211579 3.49703496
[46,] -0.01818319 -1.26211579
[47,] -6.01794552 -0.01818319
[48,] -3.21410020 -6.01794552
[49,] -0.33671835 -3.21410020
[50,] -1.67824766 -0.33671835
[51,] 0.92826945 -1.67824766
[52,] 1.32509091 0.92826945
[53,] -1.26382955 1.32509091
[54,] 2.89841712 -1.26382955
[55,] 2.64800747 2.89841712
[56,] -2.11445003 2.64800747
[57,] 3.17368008 -2.11445003
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2.28934643 -0.96985742
2 2.48133742 -2.28934643
3 -0.99720443 2.48133742
4 -0.34150497 -0.99720443
5 -0.85757508 -0.34150497
6 -1.32890930 -0.85757508
7 -2.40460004 -1.32890930
8 4.33829782 -2.40460004
9 -1.42131796 4.33829782
10 0.92312783 -1.42131796
11 0.18504947 0.92312783
12 4.26853933 0.18504947
13 1.04105413 4.26853933
14 0.96816803 1.04105413
15 4.24593307 0.96816803
16 1.83385285 4.24593307
17 0.16923357 1.83385285
18 1.27368677 0.16923357
19 1.31220692 1.27368677
20 -3.88545304 1.31220692
21 -0.84069324 -3.88545304
22 -0.85462292 -0.84069324
23 2.92614376 -0.85462292
24 -2.54572057 2.92614376
25 1.10268555 -2.54572057
26 -1.42743034 1.10268555
27 -1.56242336 -1.42743034
28 -3.39741667 -1.56242336
29 -4.23321850 -3.39741667
30 -4.38607261 -4.23321850
31 0.57567878 -4.38607261
32 -1.83542970 0.57567878
33 0.35044691 -1.83542970
34 -0.05032172 0.35044691
35 2.90675230 -0.05032172
36 2.46113885 2.90675230
37 0.48232510 2.46113885
38 -0.34382746 0.48232510
39 -2.61457472 -0.34382746
40 0.57997788 -2.61457472
41 6.18538957 0.57997788
42 1.54287801 6.18538957
43 -2.13129313 1.54287801
44 3.49703496 -2.13129313
45 -1.26211579 3.49703496
46 -0.01818319 -1.26211579
47 -6.01794552 -0.01818319
48 -3.21410020 -6.01794552
49 -0.33671835 -3.21410020
50 -1.67824766 -0.33671835
51 0.92826945 -1.67824766
52 1.32509091 0.92826945
53 -1.26382955 1.32509091
54 2.89841712 -1.26382955
55 2.64800747 2.89841712
56 -2.11445003 2.64800747
57 3.17368008 -2.11445003
> 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/7vc3c1258660722.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/8lqhd1258660722.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/9izsj1258660722.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/10f36f1258660722.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/11myd91258660722.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/128box1258660722.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/13pu2k1258660722.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/14rxk31258660722.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/15erjm1258660722.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/16md3p1258660722.tab")
+ }
>
> system("convert tmp/1vn9v1258660722.ps tmp/1vn9v1258660722.png")
> system("convert tmp/26r0t1258660722.ps tmp/26r0t1258660722.png")
> system("convert tmp/3xcq11258660722.ps tmp/3xcq11258660722.png")
> system("convert tmp/41n051258660722.ps tmp/41n051258660722.png")
> system("convert tmp/5geq91258660722.ps tmp/5geq91258660722.png")
> system("convert tmp/67wuo1258660722.ps tmp/67wuo1258660722.png")
> system("convert tmp/7vc3c1258660722.ps tmp/7vc3c1258660722.png")
> system("convert tmp/8lqhd1258660722.ps tmp/8lqhd1258660722.png")
> system("convert tmp/9izsj1258660722.ps tmp/9izsj1258660722.png")
> system("convert tmp/10f36f1258660722.ps tmp/10f36f1258660722.png")
>
>
> proc.time()
user system elapsed
2.399 1.592 2.801