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(107.11,107.56,107.57,107.70,107.81,107.67,108.75,107.67,109.43,107.72,109.62,108.35,109.54,108.25,109.53,108.26,109.84,108.31,109.67,108.33,109.79,108.36,109.56,108.36,110.22,108.97,110.40,109.62,110.69,109.60,110.72,109.64,110.89,109.65,110.58,109.64,110.94,109.93,110.91,109.81,111.22,109.77,111.09,110.10,111.00,110.40,111.06,110.50,111.55,111.89,112.32,112.10,112.64,111.92,112.36,112.15,112.04,112.16,112.37,112.17,112.59,112.32,112.89,112.38,113.22,112.34,112.85,113.14,113.06,113.18,112.99,113.21,113.32,113.76,113.74,113.99,113.91,113.95,114.52,113.93,114.96,114.01,114.91,114.10,115.30,114.11,115.44,114.10,115.52,114.12,116.08,114.68,115.94,114.71,115.56,114.73,115.88,115.81,116.66,116.01,117.41,116.12,117.68,116.49,117.85,116.51,118.21,116.60,118.92,117.01,119.03,117.01,119.17,117.12,118.95,117.22,118.92,118.38,118.90,118.80),dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> 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 = '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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 107.11 107.56 1 0 0 0 0 0 0 0 0 0 0
2 107.57 107.70 0 1 0 0 0 0 0 0 0 0 0
3 107.81 107.67 0 0 1 0 0 0 0 0 0 0 0
4 108.75 107.67 0 0 0 1 0 0 0 0 0 0 0
5 109.43 107.72 0 0 0 0 1 0 0 0 0 0 0
6 109.62 108.35 0 0 0 0 0 1 0 0 0 0 0
7 109.54 108.25 0 0 0 0 0 0 1 0 0 0 0
8 109.53 108.26 0 0 0 0 0 0 0 1 0 0 0
9 109.84 108.31 0 0 0 0 0 0 0 0 1 0 0
10 109.67 108.33 0 0 0 0 0 0 0 0 0 1 0
11 109.79 108.36 0 0 0 0 0 0 0 0 0 0 1
12 109.56 108.36 0 0 0 0 0 0 0 0 0 0 0
13 110.22 108.97 1 0 0 0 0 0 0 0 0 0 0
14 110.40 109.62 0 1 0 0 0 0 0 0 0 0 0
15 110.69 109.60 0 0 1 0 0 0 0 0 0 0 0
16 110.72 109.64 0 0 0 1 0 0 0 0 0 0 0
17 110.89 109.65 0 0 0 0 1 0 0 0 0 0 0
18 110.58 109.64 0 0 0 0 0 1 0 0 0 0 0
19 110.94 109.93 0 0 0 0 0 0 1 0 0 0 0
20 110.91 109.81 0 0 0 0 0 0 0 1 0 0 0
21 111.22 109.77 0 0 0 0 0 0 0 0 1 0 0
22 111.09 110.10 0 0 0 0 0 0 0 0 0 1 0
23 111.00 110.40 0 0 0 0 0 0 0 0 0 0 1
24 111.06 110.50 0 0 0 0 0 0 0 0 0 0 0
25 111.55 111.89 1 0 0 0 0 0 0 0 0 0 0
26 112.32 112.10 0 1 0 0 0 0 0 0 0 0 0
27 112.64 111.92 0 0 1 0 0 0 0 0 0 0 0
28 112.36 112.15 0 0 0 1 0 0 0 0 0 0 0
29 112.04 112.16 0 0 0 0 1 0 0 0 0 0 0
30 112.37 112.17 0 0 0 0 0 1 0 0 0 0 0
31 112.59 112.32 0 0 0 0 0 0 1 0 0 0 0
32 112.89 112.38 0 0 0 0 0 0 0 1 0 0 0
33 113.22 112.34 0 0 0 0 0 0 0 0 1 0 0
34 112.85 113.14 0 0 0 0 0 0 0 0 0 1 0
35 113.06 113.18 0 0 0 0 0 0 0 0 0 0 1
36 112.99 113.21 0 0 0 0 0 0 0 0 0 0 0
37 113.32 113.76 1 0 0 0 0 0 0 0 0 0 0
38 113.74 113.99 0 1 0 0 0 0 0 0 0 0 0
39 113.91 113.95 0 0 1 0 0 0 0 0 0 0 0
40 114.52 113.93 0 0 0 1 0 0 0 0 0 0 0
41 114.96 114.01 0 0 0 0 1 0 0 0 0 0 0
42 114.91 114.10 0 0 0 0 0 1 0 0 0 0 0
43 115.30 114.11 0 0 0 0 0 0 1 0 0 0 0
44 115.44 114.10 0 0 0 0 0 0 0 1 0 0 0
45 115.52 114.12 0 0 0 0 0 0 0 0 1 0 0
46 116.08 114.68 0 0 0 0 0 0 0 0 0 1 0
47 115.94 114.71 0 0 0 0 0 0 0 0 0 0 1
48 115.56 114.73 0 0 0 0 0 0 0 0 0 0 0
49 115.88 115.81 1 0 0 0 0 0 0 0 0 0 0
50 116.66 116.01 0 1 0 0 0 0 0 0 0 0 0
51 117.41 116.12 0 0 1 0 0 0 0 0 0 0 0
52 117.68 116.49 0 0 0 1 0 0 0 0 0 0 0
53 117.85 116.51 0 0 0 0 1 0 0 0 0 0 0
54 118.21 116.60 0 0 0 0 0 1 0 0 0 0 0
55 118.92 117.01 0 0 0 0 0 0 1 0 0 0 0
56 119.03 117.01 0 0 0 0 0 0 0 1 0 0 0
57 119.17 117.12 0 0 0 0 0 0 0 0 1 0 0
58 118.95 117.22 0 0 0 0 0 0 0 0 0 1 0
59 118.92 118.38 0 0 0 0 0 0 0 0 0 0 1
60 118.90 118.80 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
-0.0957 1.0052 -0.4681 -0.2336 0.1526 0.3420
M5 M6 M7 M8 M9 M10
0.5358 0.4769 0.6441 0.7582 0.9721 0.5422
M11
0.2426
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.32632 -0.37754 0.03658 0.34126 1.24570
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.09569 2.76081 -0.035 0.9725
X 1.00521 0.02429 41.378 <2e-16 ***
M1 -0.46807 0.37663 -1.243 0.2201
M2 -0.23356 0.37602 -0.621 0.5375
M3 0.15261 0.37608 0.406 0.6867
M4 0.34196 0.37584 0.910 0.3675
M5 0.53579 0.37578 1.426 0.1605
M6 0.47694 0.37552 1.270 0.2103
M7 0.64415 0.37531 1.716 0.0927 .
M8 0.75821 0.37533 2.020 0.0491 *
M9 0.97211 0.37530 2.590 0.0127 *
M10 0.54222 0.37496 1.446 0.1548
M11 0.24259 0.37482 0.647 0.5206
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5926 on 47 degrees of freedom
Multiple R-squared: 0.9744, Adjusted R-squared: 0.9679
F-statistic: 149.2 on 12 and 47 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.6271819 0.7456363 0.3728181
[2,] 0.7606695 0.4786610 0.2393305
[3,] 0.7179153 0.5641694 0.2820847
[4,] 0.6654075 0.6691850 0.3345925
[5,] 0.5812312 0.8375376 0.4187688
[6,] 0.4990498 0.9980997 0.5009502
[7,] 0.4697767 0.9395534 0.5302233
[8,] 0.5468287 0.9063426 0.4531713
[9,] 0.6489340 0.7021320 0.3510660
[10,] 0.6271838 0.7456324 0.3728162
[11,] 0.5717270 0.8565460 0.4282730
[12,] 0.5554378 0.8891244 0.4445622
[13,] 0.5219958 0.9560084 0.4780042
[14,] 0.6478995 0.7042011 0.3521005
[15,] 0.5863030 0.8273940 0.4136970
[16,] 0.5451145 0.9097711 0.4548855
[17,] 0.4828472 0.9656944 0.5171528
[18,] 0.3903021 0.7806042 0.6096979
[19,] 0.6395004 0.7209992 0.3604996
[20,] 0.6054654 0.7890692 0.3945346
[21,] 0.5221745 0.9556511 0.4778255
[22,] 0.4313711 0.8627423 0.5686289
[23,] 0.3987200 0.7974400 0.6012800
[24,] 0.5186710 0.9626581 0.4813290
[25,] 0.4615148 0.9230296 0.5384852
[26,] 0.3883342 0.7766684 0.6116658
[27,] 0.3917398 0.7834796 0.6082602
[28,] 0.4118632 0.8237264 0.5881368
[29,] 0.4532546 0.9065093 0.5467454
> postscript(file="/var/www/html/rcomp/tmp/106171258752191.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/2mzt11258752191.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/3ksog1258752191.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/4n1y31258752191.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/5yhtr1258752191.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 = 60
Frequency = 1
1 2 3 4 5 6
-0.446950020 -0.362188926 -0.478199352 0.272447057 0.708363649 0.323923978
7 8 9 10 11 12
0.177237647 0.043122962 0.088966573 0.328749409 0.718219467 0.730813746
13 14 15 16 17 18
1.245699690 0.537802168 0.461739613 0.262177502 0.228302613 -0.012800755
19 20 21 22 23 24
-0.111520145 -0.134957145 0.001355634 -0.030477551 -0.122414995 0.079657986
25 26 27 28 29 30
-0.359522188 -0.035126002 0.079645518 -0.620907057 -1.144781946 -0.765989574
31 32 33 34 35 36
-0.863979148 -0.738354482 -0.582041704 -1.326324985 -0.856907057 -0.714469167
37 38 39 40 41 42
-0.469270445 -0.514978519 -0.690936815 -0.250186147 -0.084425944 -0.166050610
43 44 45 46 47 48
0.046689632 0.082679206 -0.071320794 0.355647038 0.485117096 0.327607115
49 50 51 52 53 54
0.030042963 0.374491278 0.627751036 0.336468645 0.292541627 0.620916961
55 56 57 58 59 60
0.751572014 0.747509458 0.563040291 0.672406090 -0.224014511 -0.423609680
> postscript(file="/var/www/html/rcomp/tmp/6rhni1258752191.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.446950020 NA
1 -0.362188926 -0.446950020
2 -0.478199352 -0.362188926
3 0.272447057 -0.478199352
4 0.708363649 0.272447057
5 0.323923978 0.708363649
6 0.177237647 0.323923978
7 0.043122962 0.177237647
8 0.088966573 0.043122962
9 0.328749409 0.088966573
10 0.718219467 0.328749409
11 0.730813746 0.718219467
12 1.245699690 0.730813746
13 0.537802168 1.245699690
14 0.461739613 0.537802168
15 0.262177502 0.461739613
16 0.228302613 0.262177502
17 -0.012800755 0.228302613
18 -0.111520145 -0.012800755
19 -0.134957145 -0.111520145
20 0.001355634 -0.134957145
21 -0.030477551 0.001355634
22 -0.122414995 -0.030477551
23 0.079657986 -0.122414995
24 -0.359522188 0.079657986
25 -0.035126002 -0.359522188
26 0.079645518 -0.035126002
27 -0.620907057 0.079645518
28 -1.144781946 -0.620907057
29 -0.765989574 -1.144781946
30 -0.863979148 -0.765989574
31 -0.738354482 -0.863979148
32 -0.582041704 -0.738354482
33 -1.326324985 -0.582041704
34 -0.856907057 -1.326324985
35 -0.714469167 -0.856907057
36 -0.469270445 -0.714469167
37 -0.514978519 -0.469270445
38 -0.690936815 -0.514978519
39 -0.250186147 -0.690936815
40 -0.084425944 -0.250186147
41 -0.166050610 -0.084425944
42 0.046689632 -0.166050610
43 0.082679206 0.046689632
44 -0.071320794 0.082679206
45 0.355647038 -0.071320794
46 0.485117096 0.355647038
47 0.327607115 0.485117096
48 0.030042963 0.327607115
49 0.374491278 0.030042963
50 0.627751036 0.374491278
51 0.336468645 0.627751036
52 0.292541627 0.336468645
53 0.620916961 0.292541627
54 0.751572014 0.620916961
55 0.747509458 0.751572014
56 0.563040291 0.747509458
57 0.672406090 0.563040291
58 -0.224014511 0.672406090
59 -0.423609680 -0.224014511
60 NA -0.423609680
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.362188926 -0.446950020
[2,] -0.478199352 -0.362188926
[3,] 0.272447057 -0.478199352
[4,] 0.708363649 0.272447057
[5,] 0.323923978 0.708363649
[6,] 0.177237647 0.323923978
[7,] 0.043122962 0.177237647
[8,] 0.088966573 0.043122962
[9,] 0.328749409 0.088966573
[10,] 0.718219467 0.328749409
[11,] 0.730813746 0.718219467
[12,] 1.245699690 0.730813746
[13,] 0.537802168 1.245699690
[14,] 0.461739613 0.537802168
[15,] 0.262177502 0.461739613
[16,] 0.228302613 0.262177502
[17,] -0.012800755 0.228302613
[18,] -0.111520145 -0.012800755
[19,] -0.134957145 -0.111520145
[20,] 0.001355634 -0.134957145
[21,] -0.030477551 0.001355634
[22,] -0.122414995 -0.030477551
[23,] 0.079657986 -0.122414995
[24,] -0.359522188 0.079657986
[25,] -0.035126002 -0.359522188
[26,] 0.079645518 -0.035126002
[27,] -0.620907057 0.079645518
[28,] -1.144781946 -0.620907057
[29,] -0.765989574 -1.144781946
[30,] -0.863979148 -0.765989574
[31,] -0.738354482 -0.863979148
[32,] -0.582041704 -0.738354482
[33,] -1.326324985 -0.582041704
[34,] -0.856907057 -1.326324985
[35,] -0.714469167 -0.856907057
[36,] -0.469270445 -0.714469167
[37,] -0.514978519 -0.469270445
[38,] -0.690936815 -0.514978519
[39,] -0.250186147 -0.690936815
[40,] -0.084425944 -0.250186147
[41,] -0.166050610 -0.084425944
[42,] 0.046689632 -0.166050610
[43,] 0.082679206 0.046689632
[44,] -0.071320794 0.082679206
[45,] 0.355647038 -0.071320794
[46,] 0.485117096 0.355647038
[47,] 0.327607115 0.485117096
[48,] 0.030042963 0.327607115
[49,] 0.374491278 0.030042963
[50,] 0.627751036 0.374491278
[51,] 0.336468645 0.627751036
[52,] 0.292541627 0.336468645
[53,] 0.620916961 0.292541627
[54,] 0.751572014 0.620916961
[55,] 0.747509458 0.751572014
[56,] 0.563040291 0.747509458
[57,] 0.672406090 0.563040291
[58,] -0.224014511 0.672406090
[59,] -0.423609680 -0.224014511
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.362188926 -0.446950020
2 -0.478199352 -0.362188926
3 0.272447057 -0.478199352
4 0.708363649 0.272447057
5 0.323923978 0.708363649
6 0.177237647 0.323923978
7 0.043122962 0.177237647
8 0.088966573 0.043122962
9 0.328749409 0.088966573
10 0.718219467 0.328749409
11 0.730813746 0.718219467
12 1.245699690 0.730813746
13 0.537802168 1.245699690
14 0.461739613 0.537802168
15 0.262177502 0.461739613
16 0.228302613 0.262177502
17 -0.012800755 0.228302613
18 -0.111520145 -0.012800755
19 -0.134957145 -0.111520145
20 0.001355634 -0.134957145
21 -0.030477551 0.001355634
22 -0.122414995 -0.030477551
23 0.079657986 -0.122414995
24 -0.359522188 0.079657986
25 -0.035126002 -0.359522188
26 0.079645518 -0.035126002
27 -0.620907057 0.079645518
28 -1.144781946 -0.620907057
29 -0.765989574 -1.144781946
30 -0.863979148 -0.765989574
31 -0.738354482 -0.863979148
32 -0.582041704 -0.738354482
33 -1.326324985 -0.582041704
34 -0.856907057 -1.326324985
35 -0.714469167 -0.856907057
36 -0.469270445 -0.714469167
37 -0.514978519 -0.469270445
38 -0.690936815 -0.514978519
39 -0.250186147 -0.690936815
40 -0.084425944 -0.250186147
41 -0.166050610 -0.084425944
42 0.046689632 -0.166050610
43 0.082679206 0.046689632
44 -0.071320794 0.082679206
45 0.355647038 -0.071320794
46 0.485117096 0.355647038
47 0.327607115 0.485117096
48 0.030042963 0.327607115
49 0.374491278 0.030042963
50 0.627751036 0.374491278
51 0.336468645 0.627751036
52 0.292541627 0.336468645
53 0.620916961 0.292541627
54 0.751572014 0.620916961
55 0.747509458 0.751572014
56 0.563040291 0.747509458
57 0.672406090 0.563040291
58 -0.224014511 0.672406090
59 -0.423609680 -0.224014511
> 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/7qjwr1258752191.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/8b2xw1258752191.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/9c5mh1258752191.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/10h4071258752191.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/11sy7m1258752191.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/12minp1258752191.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/132e4u1258752191.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/14vtdt1258752191.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/15mojb1258752191.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/16d5iq1258752191.tab")
+ }
> system("convert tmp/106171258752191.ps tmp/106171258752191.png")
> system("convert tmp/2mzt11258752191.ps tmp/2mzt11258752191.png")
> system("convert tmp/3ksog1258752191.ps tmp/3ksog1258752191.png")
> system("convert tmp/4n1y31258752191.ps tmp/4n1y31258752191.png")
> system("convert tmp/5yhtr1258752191.ps tmp/5yhtr1258752191.png")
> system("convert tmp/6rhni1258752191.ps tmp/6rhni1258752191.png")
> system("convert tmp/7qjwr1258752191.ps tmp/7qjwr1258752191.png")
> system("convert tmp/8b2xw1258752191.ps tmp/8b2xw1258752191.png")
> system("convert tmp/9c5mh1258752191.ps tmp/9c5mh1258752191.png")
> system("convert tmp/10h4071258752191.ps tmp/10h4071258752191.png")
>
>
> proc.time()
user system elapsed
2.348 1.537 3.470