R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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.
Natural language support but running in an English locale
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(6654.00
+ ,5712.00
+ ,3.30
+ ,38.60
+ ,645.00
+ ,3.00
+ ,1.00
+ ,6.60
+ ,8.30
+ ,4.50
+ ,42.00
+ ,3.00
+ ,3.39
+ ,44.50
+ ,12.50
+ ,14.00
+ ,60.00
+ ,1.00
+ ,0.92
+ ,5.70
+ ,16.50
+ ,25.00
+ ,5.00
+ ,2547.00
+ ,4603.00
+ ,3.90
+ ,69.00
+ ,624.00
+ ,3.00
+ ,10.55
+ ,179.50
+ ,9.80
+ ,27.00
+ ,180.00
+ ,4.00
+ ,0.02
+ ,0.30
+ ,19.70
+ ,19.00
+ ,35.00
+ ,1.00
+ ,160.00
+ ,169.00
+ ,6.20
+ ,30.40
+ ,392.00
+ ,4.00
+ ,3.30
+ ,25.60
+ ,14.50
+ ,28.00
+ ,63.00
+ ,1.00
+ ,52.16
+ ,440.00
+ ,9.70
+ ,50.00
+ ,230.00
+ ,1.00
+ ,0.43
+ ,6.40
+ ,12.50
+ ,7.00
+ ,112.00
+ ,5.00
+ ,465.00
+ ,423.00
+ ,3.90
+ ,30.00
+ ,281.00
+ ,5.00
+ ,0.55
+ ,2.40
+ ,10.30
+ ,2.00
+ ,187.10
+ ,419.00
+ ,3.10
+ ,40.00
+ ,365.00
+ ,5.00
+ ,0.08
+ ,1.20
+ ,8.40
+ ,3.50
+ ,42.00
+ ,1.00
+ ,3.00
+ ,25.00
+ ,8.60
+ ,50.00
+ ,28.00
+ ,2.00
+ ,0.79
+ ,3.50
+ ,10.70
+ ,6.00
+ ,42.00
+ ,2.00
+ ,0.20
+ ,5.00
+ ,10.70
+ ,10.40
+ ,120.00
+ ,2.00
+ ,1.41
+ ,17.50
+ ,6.10
+ ,34.00
+ ,1.00
+ ,60.00
+ ,81.00
+ ,18.10
+ ,7.00
+ ,1.00
+ ,529.00
+ ,680.00
+ ,28.00
+ ,400.00
+ ,5.00
+ ,27.66
+ ,115.00
+ ,3.80
+ ,20.00
+ ,148.00
+ ,5.00
+ ,0.12
+ ,1.00
+ ,14.40
+ ,3.90
+ ,16.00
+ ,3.00
+ ,207.00
+ ,406.00
+ ,12.00
+ ,39.30
+ ,252.00
+ ,1.00
+ ,85.00
+ ,325.00
+ ,6.20
+ ,41.00
+ ,310.00
+ ,1.00
+ ,36.33
+ ,119.50
+ ,13.00
+ ,16.20
+ ,63.00
+ ,1.00
+ ,0.10
+ ,4.00
+ ,13.80
+ ,9.00
+ ,28.00
+ ,5.00
+ ,1.04
+ ,5.50
+ ,8.20
+ ,7.60
+ ,68.00
+ ,5.00
+ ,521.00
+ ,655.00
+ ,2.90
+ ,46.00
+ ,336.00
+ ,5.00
+ ,100.00
+ ,157.00
+ ,10.80
+ ,22.40
+ ,100.00
+ ,1.00
+ ,35.00
+ ,56.00
+ ,16.30
+ ,33.00
+ ,3.00
+ ,0.01
+ ,0.14
+ ,9.10
+ ,2.60
+ ,21.50
+ ,5.00
+ ,0.01
+ ,0.25
+ ,19.90
+ ,24.00
+ ,50.00
+ ,1.00
+ ,62.00
+ ,1320.00
+ ,8.00
+ ,100.00
+ ,267.00
+ ,1.00
+ ,0.12
+ ,3.00
+ ,10.60
+ ,30.00
+ ,2.00
+ ,1.35
+ ,8.10
+ ,11.20
+ ,45.00
+ ,3.00
+ ,0.02
+ ,0.40
+ ,13.20
+ ,3.20
+ ,19.00
+ ,4.00
+ ,0.05
+ ,0.33
+ ,12.80
+ ,2.00
+ ,30.00
+ ,4.00
+ ,1.70
+ ,6.30
+ ,19.40
+ ,5.00
+ ,12.00
+ ,2.00
+ ,3.50
+ ,10.80
+ ,17.40
+ ,6.50
+ ,120.00
+ ,2.00
+ ,250.00
+ ,490.00
+ ,23.60
+ ,440.00
+ ,5.00
+ ,0.48
+ ,15.50
+ ,17.00
+ ,12.00
+ ,140.00
+ ,2.00
+ ,10.00
+ ,115.00
+ ,10.90
+ ,20.20
+ ,170.00
+ ,4.00
+ ,1.62
+ ,11.40
+ ,13.70
+ ,13.00
+ ,17.00
+ ,2.00
+ ,192.00
+ ,180.00
+ ,8.40
+ ,27.00
+ ,115.00
+ ,4.00
+ ,2.50
+ ,12.10
+ ,8.40
+ ,18.00
+ ,31.00
+ ,5.00
+ ,4.29
+ ,39.20
+ ,12.50
+ ,13.70
+ ,63.00
+ ,2.00
+ ,0.28
+ ,1.90
+ ,13.20
+ ,4.70
+ ,21.00
+ ,3.00
+ ,4.24
+ ,50.40
+ ,9.80
+ ,9.80
+ ,52.00
+ ,1.00
+ ,6.80
+ ,179.00
+ ,9.60
+ ,29.00
+ ,164.00
+ ,2.00
+ ,0.75
+ ,12.30
+ ,6.60
+ ,7.00
+ ,225.00
+ ,2.00
+ ,3.60
+ ,21.00
+ ,5.40
+ ,6.00
+ ,225.00
+ ,3.00
+ ,14.83
+ ,98.20
+ ,2.60
+ ,17.00
+ ,150.00
+ ,5.00
+ ,55.50
+ ,175.00
+ ,3.80
+ ,20.00
+ ,151.00
+ ,5.00
+ ,1.40
+ ,12.50
+ ,11.00
+ ,12.70
+ ,90.00
+ ,2.00
+ ,0.06
+ ,1.00
+ ,10.30
+ ,3.50
+ ,3.00
+ ,0.90
+ ,2.60
+ ,13.30
+ ,4.50
+ ,60.00
+ ,2.00
+ ,2.00
+ ,12.30
+ ,5.40
+ ,7.50
+ ,200.00
+ ,3.00
+ ,0.10
+ ,2.50
+ ,15.80
+ ,2.30
+ ,46.00
+ ,3.00
+ ,4.19
+ ,58.00
+ ,10.30
+ ,24.00
+ ,210.00
+ ,4.00
+ ,3.50
+ ,3.90
+ ,19.40
+ ,3.00
+ ,14.00
+ ,2.00
+ ,4.05
+ ,17.00
+ ,13.00
+ ,38.00
+ ,3.00)
+ ,dim=c(6
+ ,62)
+ ,dimnames=list(c('G'
+ ,'H'
+ ,'Y'
+ ,'J'
+ ,'Z'
+ ,'P')
+ ,1:62))
> y <- array(NA,dim=c(6,62),dimnames=list(c('G','H','Y','J','Z','P'),1:62))
> 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.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 G H J Z P
1 3.30 6654.00 5712.00 38.60 645.00 3.00
2 8.30 1.00 6.60 4.50 42.00 3.00
3 12.50 3.39 44.50 14.00 60.00 1.00
4 16.50 0.92 5.70 25.00 5.00 2547.00
5 69.00 4603.00 3.90 624.00 3.00 10.55
6 27.00 179.50 9.80 180.00 4.00 0.02
7 19.00 0.30 19.70 35.00 1.00 160.00
8 30.40 169.00 6.20 392.00 4.00 3.30
9 28.00 25.60 14.50 63.00 1.00 52.16
10 50.00 440.00 9.70 230.00 1.00 0.43
11 7.00 6.40 12.50 112.00 5.00 465.00
12 30.00 423.00 3.90 281.00 5.00 0.55
13 2.00 2.40 10.30 187.10 419.00 3.10
14 5.00 40.00 365.00 0.08 1.20 8.40
15 1.00 3.50 42.00 3.00 25.00 8.60
16 2.00 50.00 28.00 0.79 3.50 10.70
17 2.00 6.00 42.00 0.20 5.00 10.70
18 2.00 10.40 120.00 1.41 17.50 6.10
19 60.00 34.00 1.00 81.00 18.10 7.00
20 680.00 1.00 529.00 28.00 400.00 5.00
21 3.80 27.66 115.00 20.00 148.00 5.00
22 14.40 0.12 1.00 3.90 16.00 3.00
23 12.00 207.00 406.00 39.30 252.00 1.00
24 6.20 85.00 325.00 41.00 310.00 1.00
25 13.00 36.33 119.50 16.20 63.00 1.00
26 13.80 0.10 4.00 9.00 28.00 5.00
27 8.20 1.04 5.50 7.60 68.00 5.00
28 2.90 521.00 655.00 46.00 336.00 5.00
29 10.80 100.00 157.00 22.40 100.00 1.00
30 16.30 35.00 56.00 33.00 3.00 0.01
31 2.60 0.14 9.10 21.50 5.00 0.01
32 24.00 0.25 19.90 50.00 1.00 62.00
33 100.00 1320.00 8.00 267.00 1.00 0.12
34 30.00 3.00 10.60 2.00 1.35 8.10
35 3.00 11.20 45.00 0.02 0.40 13.20
36 4.00 3.20 19.00 0.05 0.33 12.80
37 4.00 2.00 30.00 1.70 6.30 19.40
38 2.00 5.00 12.00 3.50 10.80 17.40
39 2.00 6.50 120.00 250.00 490.00 23.60
40 0.48 440.00 5.00 15.50 17.00 12.00
41 10.00 140.00 2.00 115.00 10.90 20.20
42 1.62 170.00 4.00 11.40 13.70 13.00
43 192.00 17.00 2.00 180.00 8.40 27.00
44 2.50 115.00 4.00 12.10 8.40 18.00
45 4.29 31.00 5.00 39.20 12.50 13.70
46 0.28 63.00 2.00 1.90 13.20 4.70
47 4.24 21.00 3.00 50.40 9.80 9.80
48 6.80 52.00 1.00 179.00 9.60 29.00
49 0.75 164.00 2.00 12.30 6.60 7.00
50 3.60 225.00 2.00 21.00 5.40 6.00
51 14.83 225.00 3.00 98.20 2.60 17.00
52 55.50 150.00 5.00 175.00 3.80 20.00
53 1.40 151.00 5.00 12.50 11.00 12.70
54 0.06 90.00 2.00 1.00 10.30 3.50
55 2.60 3.00 0.90 13.30 4.50 60.00
56 12.30 2.00 2.00 5.40 7.50 200.00
57 2.50 3.00 0.10 15.80 2.30 46.00
58 58.00 3.00 4.19 10.30 24.00 210.00
59 3.90 4.00 3.50 19.40 3.00 14.00
60 17.00 2.00 4.05 13.00 38.00 3.00
61 3.30 6654.00 5712.00 38.60 645.00 3.00
62 8.30 1.00 6.60 4.50 42.00 3.00
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) G H J Z P
10.1204183 -0.0120504 -0.0130610 0.1197117 0.2452467 0.0009898
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-156.5969 -13.8736 -10.1558 -0.3198 575.3454
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.1204183 14.7454741 0.686 0.4953
G -0.0120504 0.0286894 -0.420 0.6761
H -0.0130610 0.0421023 -0.310 0.7575
J 0.1197117 0.1588773 0.753 0.4543
Z 0.2452467 0.1221295 2.008 0.0495 *
P 0.0009898 0.0346052 0.029 0.9773
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 88.05 on 56 degrees of freedom
Multiple R-squared: 0.1067, Adjusted R-squared: 0.02699
F-statistic: 1.338 on 5 and 56 DF, p-value: 0.2616
> 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,] 1.603362e-04 3.206723e-04 9.998397e-01
[2,] 2.224800e-04 4.449600e-04 9.997775e-01
[3,] 7.991512e-05 1.598302e-04 9.999201e-01
[4,] 8.441969e-06 1.688394e-05 9.999916e-01
[5,] 1.407273e-06 2.814547e-06 9.999986e-01
[6,] 2.314876e-07 4.629753e-07 9.999998e-01
[7,] 4.214746e-08 8.429491e-08 1.000000e+00
[8,] 6.405042e-09 1.281008e-08 1.000000e+00
[9,] 8.381973e-10 1.676395e-09 1.000000e+00
[10,] 9.781815e-11 1.956363e-10 1.000000e+00
[11,] 4.246858e-09 8.493717e-09 1.000000e+00
[12,] 1.000000e+00 9.266547e-15 4.633274e-15
[13,] 1.000000e+00 2.472789e-14 1.236394e-14
[14,] 1.000000e+00 1.212992e-13 6.064962e-14
[15,] 1.000000e+00 2.016154e-13 1.008077e-13
[16,] 1.000000e+00 2.290014e-13 1.145007e-13
[17,] 1.000000e+00 1.059393e-12 5.296964e-13
[18,] 1.000000e+00 4.723233e-12 2.361617e-12
[19,] 1.000000e+00 1.899783e-11 9.498915e-12
[20,] 1.000000e+00 2.883364e-11 1.441682e-11
[21,] 1.000000e+00 1.061483e-10 5.307417e-11
[22,] 1.000000e+00 4.570663e-10 2.285331e-10
[23,] 1.000000e+00 1.884583e-09 9.422916e-10
[24,] 1.000000e+00 7.347913e-09 3.673956e-09
[25,] 1.000000e+00 1.708343e-08 8.541715e-09
[26,] 1.000000e+00 4.891253e-08 2.445626e-08
[27,] 9.999999e-01 1.803169e-07 9.015843e-08
[28,] 9.999997e-01 6.424231e-07 3.212116e-07
[29,] 9.999989e-01 2.192322e-06 1.096161e-06
[30,] 9.999964e-01 7.177124e-06 3.588562e-06
[31,] 9.999962e-01 7.630892e-06 3.815446e-06
[32,] 9.999891e-01 2.182767e-05 1.091383e-05
[33,] 9.999766e-01 4.686193e-05 2.343096e-05
[34,] 9.999291e-01 1.418226e-04 7.091128e-05
[35,] 1.000000e+00 2.768768e-10 1.384384e-10
[36,] 1.000000e+00 2.463690e-09 1.231845e-09
[37,] 1.000000e+00 2.081236e-08 1.040618e-08
[38,] 9.999999e-01 1.664977e-07 8.324886e-08
[39,] 9.999994e-01 1.259284e-06 6.296419e-07
[40,] 9.999999e-01 2.234623e-07 1.117312e-07
[41,] 9.999989e-01 2.264271e-06 1.132135e-06
[42,] 9.999904e-01 1.912522e-05 9.562610e-06
[43,] 9.999402e-01 1.195420e-04 5.977099e-05
[44,] 9.996195e-01 7.610633e-04 3.805317e-04
[45,] 9.995353e-01 9.293790e-04 4.646895e-04
> postscript(file="/var/www/html/freestat/rcomp/tmp/1rvmi1296221601.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/www/html/freestat/rcomp/tmp/29x431296221601.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/www/html/freestat/rcomp/tmp/3d70v1296221601.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/www/html/freestat/rcomp/tmp/4pvdb1296221601.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/www/html/freestat/rcomp/tmp/5zd941296221601.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 = 62
Frequency = 1
1 2 3 4 5
-14.84067516 -12.56419777 -13.39010742 -0.27501833 38.95219357
6 7 8 9 10
-3.35848803 4.54696853 -25.51416433 10.53874239 17.52908232
11 12 13 14 15
-17.97425046 -9.83793040 -133.11645263 -0.18332898 -15.02849486
16 17 18 19 20
-8.11571779 -8.76032133 -10.89442374 36.15981452 575.34535005
21 22 23 24 25
-43.18078128 -0.09970318 -56.83104620 -79.58695115 -12.51269902
26 27 28 29 30
-4.21623063 -19.42758225 -80.30178455 -23.27200280 2.64652475
31 32 33 34 35
-11.19992113 7.85030643 73.68219838 19.47565527 -6.51126811
36 37 38 39 40
-5.93328513 -7.47225456 -10.98831262 -156.59692685 -10.30954325
41 42 43 44 45
-14.86727049 -11.13706820 158.47565630 -9.70877983 -13.16339406
46 47 48 49 50
-12.52448210 -14.03476467 -26.49220491 -10.46604270 -7.62717470
51 52 53 54 55
-4.95005465 25.35116358 -11.04218459 -11.59897781 -10.22767798
56 57 58 59 60
-0.45395566 -10.08400567 40.64364253 -9.19850802 -3.92201569
61 62
-14.84067516 -12.56419777
> postscript(file="/var/www/html/freestat/rcomp/tmp/691w71296221601.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 = 62
Frequency = 1
lag(myerror, k = 1) myerror
0 -14.84067516 NA
1 -12.56419777 -14.84067516
2 -13.39010742 -12.56419777
3 -0.27501833 -13.39010742
4 38.95219357 -0.27501833
5 -3.35848803 38.95219357
6 4.54696853 -3.35848803
7 -25.51416433 4.54696853
8 10.53874239 -25.51416433
9 17.52908232 10.53874239
10 -17.97425046 17.52908232
11 -9.83793040 -17.97425046
12 -133.11645263 -9.83793040
13 -0.18332898 -133.11645263
14 -15.02849486 -0.18332898
15 -8.11571779 -15.02849486
16 -8.76032133 -8.11571779
17 -10.89442374 -8.76032133
18 36.15981452 -10.89442374
19 575.34535005 36.15981452
20 -43.18078128 575.34535005
21 -0.09970318 -43.18078128
22 -56.83104620 -0.09970318
23 -79.58695115 -56.83104620
24 -12.51269902 -79.58695115
25 -4.21623063 -12.51269902
26 -19.42758225 -4.21623063
27 -80.30178455 -19.42758225
28 -23.27200280 -80.30178455
29 2.64652475 -23.27200280
30 -11.19992113 2.64652475
31 7.85030643 -11.19992113
32 73.68219838 7.85030643
33 19.47565527 73.68219838
34 -6.51126811 19.47565527
35 -5.93328513 -6.51126811
36 -7.47225456 -5.93328513
37 -10.98831262 -7.47225456
38 -156.59692685 -10.98831262
39 -10.30954325 -156.59692685
40 -14.86727049 -10.30954325
41 -11.13706820 -14.86727049
42 158.47565630 -11.13706820
43 -9.70877983 158.47565630
44 -13.16339406 -9.70877983
45 -12.52448210 -13.16339406
46 -14.03476467 -12.52448210
47 -26.49220491 -14.03476467
48 -10.46604270 -26.49220491
49 -7.62717470 -10.46604270
50 -4.95005465 -7.62717470
51 25.35116358 -4.95005465
52 -11.04218459 25.35116358
53 -11.59897781 -11.04218459
54 -10.22767798 -11.59897781
55 -0.45395566 -10.22767798
56 -10.08400567 -0.45395566
57 40.64364253 -10.08400567
58 -9.19850802 40.64364253
59 -3.92201569 -9.19850802
60 -14.84067516 -3.92201569
61 -12.56419777 -14.84067516
62 NA -12.56419777
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -12.56419777 -14.84067516
[2,] -13.39010742 -12.56419777
[3,] -0.27501833 -13.39010742
[4,] 38.95219357 -0.27501833
[5,] -3.35848803 38.95219357
[6,] 4.54696853 -3.35848803
[7,] -25.51416433 4.54696853
[8,] 10.53874239 -25.51416433
[9,] 17.52908232 10.53874239
[10,] -17.97425046 17.52908232
[11,] -9.83793040 -17.97425046
[12,] -133.11645263 -9.83793040
[13,] -0.18332898 -133.11645263
[14,] -15.02849486 -0.18332898
[15,] -8.11571779 -15.02849486
[16,] -8.76032133 -8.11571779
[17,] -10.89442374 -8.76032133
[18,] 36.15981452 -10.89442374
[19,] 575.34535005 36.15981452
[20,] -43.18078128 575.34535005
[21,] -0.09970318 -43.18078128
[22,] -56.83104620 -0.09970318
[23,] -79.58695115 -56.83104620
[24,] -12.51269902 -79.58695115
[25,] -4.21623063 -12.51269902
[26,] -19.42758225 -4.21623063
[27,] -80.30178455 -19.42758225
[28,] -23.27200280 -80.30178455
[29,] 2.64652475 -23.27200280
[30,] -11.19992113 2.64652475
[31,] 7.85030643 -11.19992113
[32,] 73.68219838 7.85030643
[33,] 19.47565527 73.68219838
[34,] -6.51126811 19.47565527
[35,] -5.93328513 -6.51126811
[36,] -7.47225456 -5.93328513
[37,] -10.98831262 -7.47225456
[38,] -156.59692685 -10.98831262
[39,] -10.30954325 -156.59692685
[40,] -14.86727049 -10.30954325
[41,] -11.13706820 -14.86727049
[42,] 158.47565630 -11.13706820
[43,] -9.70877983 158.47565630
[44,] -13.16339406 -9.70877983
[45,] -12.52448210 -13.16339406
[46,] -14.03476467 -12.52448210
[47,] -26.49220491 -14.03476467
[48,] -10.46604270 -26.49220491
[49,] -7.62717470 -10.46604270
[50,] -4.95005465 -7.62717470
[51,] 25.35116358 -4.95005465
[52,] -11.04218459 25.35116358
[53,] -11.59897781 -11.04218459
[54,] -10.22767798 -11.59897781
[55,] -0.45395566 -10.22767798
[56,] -10.08400567 -0.45395566
[57,] 40.64364253 -10.08400567
[58,] -9.19850802 40.64364253
[59,] -3.92201569 -9.19850802
[60,] -14.84067516 -3.92201569
[61,] -12.56419777 -14.84067516
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -12.56419777 -14.84067516
2 -13.39010742 -12.56419777
3 -0.27501833 -13.39010742
4 38.95219357 -0.27501833
5 -3.35848803 38.95219357
6 4.54696853 -3.35848803
7 -25.51416433 4.54696853
8 10.53874239 -25.51416433
9 17.52908232 10.53874239
10 -17.97425046 17.52908232
11 -9.83793040 -17.97425046
12 -133.11645263 -9.83793040
13 -0.18332898 -133.11645263
14 -15.02849486 -0.18332898
15 -8.11571779 -15.02849486
16 -8.76032133 -8.11571779
17 -10.89442374 -8.76032133
18 36.15981452 -10.89442374
19 575.34535005 36.15981452
20 -43.18078128 575.34535005
21 -0.09970318 -43.18078128
22 -56.83104620 -0.09970318
23 -79.58695115 -56.83104620
24 -12.51269902 -79.58695115
25 -4.21623063 -12.51269902
26 -19.42758225 -4.21623063
27 -80.30178455 -19.42758225
28 -23.27200280 -80.30178455
29 2.64652475 -23.27200280
30 -11.19992113 2.64652475
31 7.85030643 -11.19992113
32 73.68219838 7.85030643
33 19.47565527 73.68219838
34 -6.51126811 19.47565527
35 -5.93328513 -6.51126811
36 -7.47225456 -5.93328513
37 -10.98831262 -7.47225456
38 -156.59692685 -10.98831262
39 -10.30954325 -156.59692685
40 -14.86727049 -10.30954325
41 -11.13706820 -14.86727049
42 158.47565630 -11.13706820
43 -9.70877983 158.47565630
44 -13.16339406 -9.70877983
45 -12.52448210 -13.16339406
46 -14.03476467 -12.52448210
47 -26.49220491 -14.03476467
48 -10.46604270 -26.49220491
49 -7.62717470 -10.46604270
50 -4.95005465 -7.62717470
51 25.35116358 -4.95005465
52 -11.04218459 25.35116358
53 -11.59897781 -11.04218459
54 -10.22767798 -11.59897781
55 -0.45395566 -10.22767798
56 -10.08400567 -0.45395566
57 40.64364253 -10.08400567
58 -9.19850802 40.64364253
59 -3.92201569 -9.19850802
60 -14.84067516 -3.92201569
61 -12.56419777 -14.84067516
> 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/freestat/rcomp/tmp/7rqgg1296221601.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/www/html/freestat/rcomp/tmp/8ufzq1296221601.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/www/html/freestat/rcomp/tmp/9l62i1296221601.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/www/html/freestat/rcomp/tmp/10q5r11296221601.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/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11u1dr1296221601.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/freestat/rcomp/tmp/124pgy1296221601.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/freestat/rcomp/tmp/13mnct1296221601.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/freestat/rcomp/tmp/149zex1296221601.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/freestat/rcomp/tmp/15tcxi1296221601.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/freestat/rcomp/tmp/16fqyy1296221601.tab")
+ }
>
> try(system("convert tmp/1rvmi1296221601.ps tmp/1rvmi1296221601.png",intern=TRUE))
character(0)
> try(system("convert tmp/29x431296221601.ps tmp/29x431296221601.png",intern=TRUE))
character(0)
> try(system("convert tmp/3d70v1296221601.ps tmp/3d70v1296221601.png",intern=TRUE))
character(0)
> try(system("convert tmp/4pvdb1296221601.ps tmp/4pvdb1296221601.png",intern=TRUE))
character(0)
> try(system("convert tmp/5zd941296221601.ps tmp/5zd941296221601.png",intern=TRUE))
character(0)
> try(system("convert tmp/691w71296221601.ps tmp/691w71296221601.png",intern=TRUE))
character(0)
> try(system("convert tmp/7rqgg1296221601.ps tmp/7rqgg1296221601.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ufzq1296221601.ps tmp/8ufzq1296221601.png",intern=TRUE))
character(0)
> try(system("convert tmp/9l62i1296221601.ps tmp/9l62i1296221601.png",intern=TRUE))
character(0)
> try(system("convert tmp/10q5r11296221601.ps tmp/10q5r11296221601.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.903 2.448 4.411