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(1.4,1.9,1,1.6,-0.8,0,-2.9,-1.3,-0.7,-0.4,-0.7,-0.3,1.5,1.4,3,2.6,3.2,2.8,3.1,2.6,3.9,3.4,1,1.7,1.3,1.2,0.8,0,1.2,0,2.9,1.6,3.9,2.5,4.5,3.2,4.5,3.4,3.3,2.3,2,1.9,1.5,1.7,1,1.9,2.1,3.3,3,3.8,4,4.4,5.1,4.5,4.5,3.5,4.2,3,3.3,2.8,2.7,2.9,1.8,2.6,1.4,2.1,0.5,1.5,-0.4,1.1,0.8,1.5,0.7,1.7,1.9,2.3,2,2.3,1.1,1.9,0.9,2,0.4,1.6,0.7,1.2,2.1,1.9,2.8,2.1,3.9,2.4,3.5,2.9,2,2.5,2,2.3,1.5,2.5,2.5,2.6,3.1,2.4,2.7,2.5,2.8,2.1,2.5,2.2,3,2.7,3.2,3,2.8,3.2,2.4,3,2,2.7,1.8,2.5,1.1,1.6,-1.5,0.1,-3.7,-1.9),dim=c(2,64),dimnames=list(c('bbp','dnst'),1:64))
> y <- array(NA,dim=c(2,64),dimnames=list(c('bbp','dnst'),1:64))
> 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
bbp dnst M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 1.4 1.9 1 0 0 0 0 0 0 0 0 0 0
2 1.0 1.6 0 1 0 0 0 0 0 0 0 0 0
3 -0.8 0.0 0 0 1 0 0 0 0 0 0 0 0
4 -2.9 -1.3 0 0 0 1 0 0 0 0 0 0 0
5 -0.7 -0.4 0 0 0 0 1 0 0 0 0 0 0
6 -0.7 -0.3 0 0 0 0 0 1 0 0 0 0 0
7 1.5 1.4 0 0 0 0 0 0 1 0 0 0 0
8 3.0 2.6 0 0 0 0 0 0 0 1 0 0 0
9 3.2 2.8 0 0 0 0 0 0 0 0 1 0 0
10 3.1 2.6 0 0 0 0 0 0 0 0 0 1 0
11 3.9 3.4 0 0 0 0 0 0 0 0 0 0 1
12 1.0 1.7 0 0 0 0 0 0 0 0 0 0 0
13 1.3 1.2 1 0 0 0 0 0 0 0 0 0 0
14 0.8 0.0 0 1 0 0 0 0 0 0 0 0 0
15 1.2 0.0 0 0 1 0 0 0 0 0 0 0 0
16 2.9 1.6 0 0 0 1 0 0 0 0 0 0 0
17 3.9 2.5 0 0 0 0 1 0 0 0 0 0 0
18 4.5 3.2 0 0 0 0 0 1 0 0 0 0 0
19 4.5 3.4 0 0 0 0 0 0 1 0 0 0 0
20 3.3 2.3 0 0 0 0 0 0 0 1 0 0 0
21 2.0 1.9 0 0 0 0 0 0 0 0 1 0 0
22 1.5 1.7 0 0 0 0 0 0 0 0 0 1 0
23 1.0 1.9 0 0 0 0 0 0 0 0 0 0 1
24 2.1 3.3 0 0 0 0 0 0 0 0 0 0 0
25 3.0 3.8 1 0 0 0 0 0 0 0 0 0 0
26 4.0 4.4 0 1 0 0 0 0 0 0 0 0 0
27 5.1 4.5 0 0 1 0 0 0 0 0 0 0 0
28 4.5 3.5 0 0 0 1 0 0 0 0 0 0 0
29 4.2 3.0 0 0 0 0 1 0 0 0 0 0 0
30 3.3 2.8 0 0 0 0 0 1 0 0 0 0 0
31 2.7 2.9 0 0 0 0 0 0 1 0 0 0 0
32 1.8 2.6 0 0 0 0 0 0 0 1 0 0 0
33 1.4 2.1 0 0 0 0 0 0 0 0 1 0 0
34 0.5 1.5 0 0 0 0 0 0 0 0 0 1 0
35 -0.4 1.1 0 0 0 0 0 0 0 0 0 0 1
36 0.8 1.5 0 0 0 0 0 0 0 0 0 0 0
37 0.7 1.7 1 0 0 0 0 0 0 0 0 0 0
38 1.9 2.3 0 1 0 0 0 0 0 0 0 0 0
39 2.0 2.3 0 0 1 0 0 0 0 0 0 0 0
40 1.1 1.9 0 0 0 1 0 0 0 0 0 0 0
41 0.9 2.0 0 0 0 0 1 0 0 0 0 0 0
42 0.4 1.6 0 0 0 0 0 1 0 0 0 0 0
43 0.7 1.2 0 0 0 0 0 0 1 0 0 0 0
44 2.1 1.9 0 0 0 0 0 0 0 1 0 0 0
45 2.8 2.1 0 0 0 0 0 0 0 0 1 0 0
46 3.9 2.4 0 0 0 0 0 0 0 0 0 1 0
47 3.5 2.9 0 0 0 0 0 0 0 0 0 0 1
48 2.0 2.5 0 0 0 0 0 0 0 0 0 0 0
49 2.0 2.3 1 0 0 0 0 0 0 0 0 0 0
50 1.5 2.5 0 1 0 0 0 0 0 0 0 0 0
51 2.5 2.6 0 0 1 0 0 0 0 0 0 0 0
52 3.1 2.4 0 0 0 1 0 0 0 0 0 0 0
53 2.7 2.5 0 0 0 0 1 0 0 0 0 0 0
54 2.8 2.1 0 0 0 0 0 1 0 0 0 0 0
55 2.5 2.2 0 0 0 0 0 0 1 0 0 0 0
56 3.0 2.7 0 0 0 0 0 0 0 1 0 0 0
57 3.2 3.0 0 0 0 0 0 0 0 0 1 0 0
58 2.8 3.2 0 0 0 0 0 0 0 0 0 1 0
59 2.4 3.0 0 0 0 0 0 0 0 0 0 0 1
60 2.0 2.7 0 0 0 0 0 0 0 0 0 0 0
61 1.8 2.5 1 0 0 0 0 0 0 0 0 0 0
62 1.1 1.6 0 1 0 0 0 0 0 0 0 0 0
63 -1.5 0.1 0 0 1 0 0 0 0 0 0 0 0
64 -3.7 -1.9 0 0 0 1 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) dnst M1 M2 M3 M4
-1.4320 1.2872 0.2573 0.4885 0.8106 0.9353
M5 M6 M7 M8 M9 M10
1.1606 1.0721 0.9545 0.9570 0.8885 0.8572
M11
0.3455
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.402976 -0.420137 -0.008195 0.360278 1.821390
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.43203 0.39531 -3.623 0.000672 ***
dnst 1.28719 0.08452 15.230 < 2e-16 ***
M1 0.25730 0.46353 0.555 0.581263
M2 0.48850 0.46402 1.053 0.297417
M3 0.81064 0.46784 1.733 0.089185 .
M4 0.93527 0.47642 1.963 0.055099 .
M5 1.16062 0.48535 2.391 0.020516 *
M6 1.07211 0.48561 2.208 0.031786 *
M7 0.95446 0.48416 1.971 0.054118 .
M8 0.95702 0.48410 1.977 0.053470 .
M9 0.88851 0.48407 1.836 0.072262 .
M10 0.85723 0.48408 1.771 0.082562 .
M11 0.34554 0.48416 0.714 0.478677
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7654 on 51 degrees of freedom
Multiple R-squared: 0.8363, Adjusted R-squared: 0.7978
F-statistic: 21.71 on 12 and 51 DF, p-value: 6.006e-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.9874716 0.02505678 0.01252839
[2,] 0.9783426 0.04331474 0.02165737
[3,] 0.9603683 0.07926345 0.03963172
[4,] 0.9314306 0.13713885 0.06856943
[5,] 0.9207736 0.15845285 0.07922643
[6,] 0.8725473 0.25490541 0.12745271
[7,] 0.8090942 0.38181156 0.19090578
[8,] 0.7491846 0.50163087 0.25081544
[9,] 0.7851496 0.42970071 0.21485035
[10,] 0.8546318 0.29073634 0.14536817
[11,] 0.9018590 0.19628197 0.09814099
[12,] 0.8654338 0.26913246 0.13456623
[13,] 0.8157339 0.36853211 0.18426606
[14,] 0.8119964 0.37600730 0.18800365
[15,] 0.7486073 0.50278536 0.25139268
[16,] 0.7279252 0.54414956 0.27207478
[17,] 0.8053419 0.38931612 0.19465806
[18,] 0.8020704 0.39585925 0.19792963
[19,] 0.8055102 0.38897952 0.19448976
[20,] 0.7853955 0.42920902 0.21460451
[21,] 0.7241190 0.55176206 0.27588103
[22,] 0.6476105 0.70477896 0.35238948
[23,] 0.5628880 0.87422394 0.43711197
[24,] 0.4900358 0.98007162 0.50996419
[25,] 0.5117703 0.97645934 0.48822967
[26,] 0.5864179 0.82716422 0.41358211
[27,] 0.7154253 0.56914934 0.28457467
[28,] 0.6239046 0.75219077 0.37609538
[29,] 0.5090411 0.98191772 0.49095886
[30,] 0.4553936 0.91078710 0.54460645
[31,] 0.8464673 0.30706540 0.15353270
[32,] 0.9185989 0.16280216 0.08140108
[33,] 0.8284720 0.34305592 0.17152796
> postscript(file="/var/www/html/rcomp/tmp/1tusp1258643689.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/2ra2l1258643689.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/3vdxn1258643689.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/48gd31258643689.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/59b901258643689.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 = 64
Frequency = 1
1 2 3 4 5
0.1290646515 -0.1159761546 -0.1786095720 -0.7298807728 0.0862899744
6 7 8 9 10
0.0460828208 0.1754990427 0.1283050882 0.1393785391 0.3280979346
11 12 13 14 15
0.6100376828 0.2438041309 0.9301004196 1.7435341726 1.8213904280
16 17 18 19 20
1.3372567591 0.9534275064 0.7409039801 0.6011111337 0.8144632745
21 22 23 24 25
0.0978530982 -0.1134275064 -0.3591713855 -0.7157061963 -0.7166038620
26 27 28 29 30
-0.7201192272 -0.0709823673 0.4915882456 0.6098305291 0.0557815619
31 32 33 34 35
-0.5552918891 -1.0716949118 -0.7595856927 -0.8559887155 -0.7294162219
36 37 38 39 40
0.3012429218 -0.3134965576 -0.1170119227 -0.3391556674 -0.8489014272
41 42 43 44 45
-1.4029755164 -1.2995856927 -0.3670621664 0.1293408563 0.6404143073
46 47 48 49 50
1.3855367255 0.8536346600 0.2140489673 0.2141870697 -0.7744507136
51 52 53 54 55
-0.2253138537 0.5075015955 -0.2465724936 0.4568173300 0.1457438791
56 57 58 59 60
-0.0004143073 -0.1180602518 -0.7442184381 -0.3750847354 -0.0433898236
61 62 63 64
-0.2432517212 -0.0159761546 -1.0073289675 -0.7575644001
> postscript(file="/var/www/html/rcomp/tmp/6fxdb1258643689.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 = 64
Frequency = 1
lag(myerror, k = 1) myerror
0 0.1290646515 NA
1 -0.1159761546 0.1290646515
2 -0.1786095720 -0.1159761546
3 -0.7298807728 -0.1786095720
4 0.0862899744 -0.7298807728
5 0.0460828208 0.0862899744
6 0.1754990427 0.0460828208
7 0.1283050882 0.1754990427
8 0.1393785391 0.1283050882
9 0.3280979346 0.1393785391
10 0.6100376828 0.3280979346
11 0.2438041309 0.6100376828
12 0.9301004196 0.2438041309
13 1.7435341726 0.9301004196
14 1.8213904280 1.7435341726
15 1.3372567591 1.8213904280
16 0.9534275064 1.3372567591
17 0.7409039801 0.9534275064
18 0.6011111337 0.7409039801
19 0.8144632745 0.6011111337
20 0.0978530982 0.8144632745
21 -0.1134275064 0.0978530982
22 -0.3591713855 -0.1134275064
23 -0.7157061963 -0.3591713855
24 -0.7166038620 -0.7157061963
25 -0.7201192272 -0.7166038620
26 -0.0709823673 -0.7201192272
27 0.4915882456 -0.0709823673
28 0.6098305291 0.4915882456
29 0.0557815619 0.6098305291
30 -0.5552918891 0.0557815619
31 -1.0716949118 -0.5552918891
32 -0.7595856927 -1.0716949118
33 -0.8559887155 -0.7595856927
34 -0.7294162219 -0.8559887155
35 0.3012429218 -0.7294162219
36 -0.3134965576 0.3012429218
37 -0.1170119227 -0.3134965576
38 -0.3391556674 -0.1170119227
39 -0.8489014272 -0.3391556674
40 -1.4029755164 -0.8489014272
41 -1.2995856927 -1.4029755164
42 -0.3670621664 -1.2995856927
43 0.1293408563 -0.3670621664
44 0.6404143073 0.1293408563
45 1.3855367255 0.6404143073
46 0.8536346600 1.3855367255
47 0.2140489673 0.8536346600
48 0.2141870697 0.2140489673
49 -0.7744507136 0.2141870697
50 -0.2253138537 -0.7744507136
51 0.5075015955 -0.2253138537
52 -0.2465724936 0.5075015955
53 0.4568173300 -0.2465724936
54 0.1457438791 0.4568173300
55 -0.0004143073 0.1457438791
56 -0.1180602518 -0.0004143073
57 -0.7442184381 -0.1180602518
58 -0.3750847354 -0.7442184381
59 -0.0433898236 -0.3750847354
60 -0.2432517212 -0.0433898236
61 -0.0159761546 -0.2432517212
62 -1.0073289675 -0.0159761546
63 -0.7575644001 -1.0073289675
64 NA -0.7575644001
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.1159761546 0.1290646515
[2,] -0.1786095720 -0.1159761546
[3,] -0.7298807728 -0.1786095720
[4,] 0.0862899744 -0.7298807728
[5,] 0.0460828208 0.0862899744
[6,] 0.1754990427 0.0460828208
[7,] 0.1283050882 0.1754990427
[8,] 0.1393785391 0.1283050882
[9,] 0.3280979346 0.1393785391
[10,] 0.6100376828 0.3280979346
[11,] 0.2438041309 0.6100376828
[12,] 0.9301004196 0.2438041309
[13,] 1.7435341726 0.9301004196
[14,] 1.8213904280 1.7435341726
[15,] 1.3372567591 1.8213904280
[16,] 0.9534275064 1.3372567591
[17,] 0.7409039801 0.9534275064
[18,] 0.6011111337 0.7409039801
[19,] 0.8144632745 0.6011111337
[20,] 0.0978530982 0.8144632745
[21,] -0.1134275064 0.0978530982
[22,] -0.3591713855 -0.1134275064
[23,] -0.7157061963 -0.3591713855
[24,] -0.7166038620 -0.7157061963
[25,] -0.7201192272 -0.7166038620
[26,] -0.0709823673 -0.7201192272
[27,] 0.4915882456 -0.0709823673
[28,] 0.6098305291 0.4915882456
[29,] 0.0557815619 0.6098305291
[30,] -0.5552918891 0.0557815619
[31,] -1.0716949118 -0.5552918891
[32,] -0.7595856927 -1.0716949118
[33,] -0.8559887155 -0.7595856927
[34,] -0.7294162219 -0.8559887155
[35,] 0.3012429218 -0.7294162219
[36,] -0.3134965576 0.3012429218
[37,] -0.1170119227 -0.3134965576
[38,] -0.3391556674 -0.1170119227
[39,] -0.8489014272 -0.3391556674
[40,] -1.4029755164 -0.8489014272
[41,] -1.2995856927 -1.4029755164
[42,] -0.3670621664 -1.2995856927
[43,] 0.1293408563 -0.3670621664
[44,] 0.6404143073 0.1293408563
[45,] 1.3855367255 0.6404143073
[46,] 0.8536346600 1.3855367255
[47,] 0.2140489673 0.8536346600
[48,] 0.2141870697 0.2140489673
[49,] -0.7744507136 0.2141870697
[50,] -0.2253138537 -0.7744507136
[51,] 0.5075015955 -0.2253138537
[52,] -0.2465724936 0.5075015955
[53,] 0.4568173300 -0.2465724936
[54,] 0.1457438791 0.4568173300
[55,] -0.0004143073 0.1457438791
[56,] -0.1180602518 -0.0004143073
[57,] -0.7442184381 -0.1180602518
[58,] -0.3750847354 -0.7442184381
[59,] -0.0433898236 -0.3750847354
[60,] -0.2432517212 -0.0433898236
[61,] -0.0159761546 -0.2432517212
[62,] -1.0073289675 -0.0159761546
[63,] -0.7575644001 -1.0073289675
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.1159761546 0.1290646515
2 -0.1786095720 -0.1159761546
3 -0.7298807728 -0.1786095720
4 0.0862899744 -0.7298807728
5 0.0460828208 0.0862899744
6 0.1754990427 0.0460828208
7 0.1283050882 0.1754990427
8 0.1393785391 0.1283050882
9 0.3280979346 0.1393785391
10 0.6100376828 0.3280979346
11 0.2438041309 0.6100376828
12 0.9301004196 0.2438041309
13 1.7435341726 0.9301004196
14 1.8213904280 1.7435341726
15 1.3372567591 1.8213904280
16 0.9534275064 1.3372567591
17 0.7409039801 0.9534275064
18 0.6011111337 0.7409039801
19 0.8144632745 0.6011111337
20 0.0978530982 0.8144632745
21 -0.1134275064 0.0978530982
22 -0.3591713855 -0.1134275064
23 -0.7157061963 -0.3591713855
24 -0.7166038620 -0.7157061963
25 -0.7201192272 -0.7166038620
26 -0.0709823673 -0.7201192272
27 0.4915882456 -0.0709823673
28 0.6098305291 0.4915882456
29 0.0557815619 0.6098305291
30 -0.5552918891 0.0557815619
31 -1.0716949118 -0.5552918891
32 -0.7595856927 -1.0716949118
33 -0.8559887155 -0.7595856927
34 -0.7294162219 -0.8559887155
35 0.3012429218 -0.7294162219
36 -0.3134965576 0.3012429218
37 -0.1170119227 -0.3134965576
38 -0.3391556674 -0.1170119227
39 -0.8489014272 -0.3391556674
40 -1.4029755164 -0.8489014272
41 -1.2995856927 -1.4029755164
42 -0.3670621664 -1.2995856927
43 0.1293408563 -0.3670621664
44 0.6404143073 0.1293408563
45 1.3855367255 0.6404143073
46 0.8536346600 1.3855367255
47 0.2140489673 0.8536346600
48 0.2141870697 0.2140489673
49 -0.7744507136 0.2141870697
50 -0.2253138537 -0.7744507136
51 0.5075015955 -0.2253138537
52 -0.2465724936 0.5075015955
53 0.4568173300 -0.2465724936
54 0.1457438791 0.4568173300
55 -0.0004143073 0.1457438791
56 -0.1180602518 -0.0004143073
57 -0.7442184381 -0.1180602518
58 -0.3750847354 -0.7442184381
59 -0.0433898236 -0.3750847354
60 -0.2432517212 -0.0433898236
61 -0.0159761546 -0.2432517212
62 -1.0073289675 -0.0159761546
63 -0.7575644001 -1.0073289675
> 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/7f5wi1258643689.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/84rjo1258643690.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/9khxk1258643690.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/107qi51258643690.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/11pfqq1258643690.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/12jk7g1258643690.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/13h6dp1258643690.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/14y0ef1258643690.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/15sr4i1258643690.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/16b6m51258643690.tab")
+ }
>
> system("convert tmp/1tusp1258643689.ps tmp/1tusp1258643689.png")
> system("convert tmp/2ra2l1258643689.ps tmp/2ra2l1258643689.png")
> system("convert tmp/3vdxn1258643689.ps tmp/3vdxn1258643689.png")
> system("convert tmp/48gd31258643689.ps tmp/48gd31258643689.png")
> system("convert tmp/59b901258643689.ps tmp/59b901258643689.png")
> system("convert tmp/6fxdb1258643689.ps tmp/6fxdb1258643689.png")
> system("convert tmp/7f5wi1258643689.ps tmp/7f5wi1258643689.png")
> system("convert tmp/84rjo1258643690.ps tmp/84rjo1258643690.png")
> system("convert tmp/9khxk1258643690.ps tmp/9khxk1258643690.png")
> system("convert tmp/107qi51258643690.ps tmp/107qi51258643690.png")
>
>
> proc.time()
user system elapsed
2.512 1.613 3.143