R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-pc-linux-gnu (32-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,100
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,50
+ ,100
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,100
+ ,50
+ ,50
+ ,100
+ ,100
+ ,100
+ ,100
+ ,100
+ ,100
+ ,100
+ ,100
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,100
+ ,100
+ ,50
+ ,100
+ ,100
+ ,100
+ ,50
+ ,50
+ ,100
+ ,100
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,100
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,50
+ ,50
+ ,100
+ ,100
+ ,50
+ ,100
+ ,100
+ ,50
+ ,50
+ ,50
+ ,100
+ ,100
+ ,100
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,100
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,100
+ ,100
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,100
+ ,100
+ ,50
+ ,100
+ ,100
+ ,50
+ ,100
+ ,50
+ ,100
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,100
+ ,100
+ ,50
+ ,100
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,100
+ ,100
+ ,100
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,100
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,100
+ ,100
+ ,100
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,100
+ ,100
+ ,100
+ ,100
+ ,100
+ ,100
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,100
+ ,100
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,100
+ ,50
+ ,100
+ ,100
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,100
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,100
+ ,100
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,100
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,100
+ ,100
+ ,50
+ ,50
+ ,100
+ ,100
+ ,50
+ ,50
+ ,50
+ ,100
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,100
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,50
+ ,100
+ ,50
+ ,100
+ ,50
+ ,50
+ ,100
+ ,50
+ ,50
+ ,50)
+ ,dim=c(6
+ ,86)
+ ,dimnames=list(c('Used'
+ ,'Correctanal'
+ ,'uselimit'
+ ,'useful'
+ ,'T40'
+ ,'outcome')
+ ,1:86))
> y <- array(NA,dim=c(6,86),dimnames=list(c('Used','Correctanal','uselimit','useful','T40','outcome'),1:86))
> 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 = '1'
> par3 <- 'No Linear Trend'
> par2 <- 'Do not include Seasonal 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, 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
Used Correctanal uselimit useful T40 outcome
1 50 50 100 50 100 100
2 50 50 50 50 50 50
3 50 50 50 50 50 50
4 50 50 50 50 50 50
5 50 50 50 50 50 50
6 50 50 100 100 50 100
7 50 50 50 50 50 50
8 50 50 50 50 100 50
9 50 50 50 50 50 100
10 50 50 100 50 50 50
11 50 50 100 50 100 50
12 50 50 50 50 50 50
13 100 50 50 100 50 50
14 50 50 100 50 100 50
15 100 50 50 100 50 100
16 100 50 50 100 100 100
17 100 100 100 100 100 50
18 50 50 100 50 100 50
19 50 50 50 50 50 100
20 100 100 50 100 100 100
21 50 50 100 100 50 50
22 100 50 100 100 50 100
23 50 50 50 100 50 100
24 50 50 100 100 50 100
25 100 50 50 50 100 100
26 100 50 50 100 50 50
27 50 50 100 50 50 100
28 100 50 50 50 50 50
29 50 50 50 50 50 100
30 50 50 50 100 50 50
31 50 50 50 50 50 50
32 50 50 100 50 50 50
33 50 50 100 100 50 50
34 50 50 50 50 100 100
35 50 50 50 50 50 50
36 50 50 50 50 50 50
37 100 50 100 100 100 50
38 100 50 50 50 50 100
39 50 50 50 100 50 100
40 50 50 50 100 100 50
41 100 100 50 100 50 100
42 100 50 50 50 50 100
43 50 50 100 100 50 100
44 50 50 100 50 100 50
45 50 50 50 100 50 50
46 50 50 50 100 50 100
47 50 50 50 50 50 50
48 50 50 50 50 50 100
49 50 50 50 100 50 100
50 50 50 50 50 50 50
51 100 50 50 50 100 50
52 100 100 100 100 100 50
53 50 50 50 50 50 100
54 100 100 50 50 50 50
55 50 50 50 50 50 50
56 100 50 50 50 100 100
57 100 50 50 100 50 100
58 50 50 50 50 50 100
59 50 50 50 50 50 100
60 100 100 100 100 100 100
61 50 50 100 50 100 100
62 100 50 50 100 50 50
63 50 50 50 50 50 50
64 50 50 100 50 100 100
65 50 50 50 50 50 50
66 50 50 50 50 50 50
67 100 100 50 100 100 50
68 50 50 100 50 50 50
69 50 50 50 50 50 100
70 100 50 50 50 50 50
71 50 50 50 50 50 50
72 50 50 50 50 50 100
73 100 50 50 50 50 100
74 100 50 100 50 50 50
75 50 50 50 50 50 100
76 50 50 50 100 100 100
77 50 50 50 50 50 100
78 100 50 50 100 50 100
79 100 100 50 50 100 100
80 50 50 50 100 100 50
81 50 50 50 50 50 50
82 100 50 100 50 50 100
83 50 50 50 50 50 50
84 100 100 50 50 50 50
85 50 50 50 100 50 100
86 50 50 100 50 50 50
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Correctanal uselimit useful T40 outcome
14.89985 0.67242 -0.09514 0.15932 0.07096 0.06938
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-23.730 -12.216 -8.747 0.875 46.010
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.89985 12.94798 1.151 0.253
Correctanal 0.67242 0.15600 4.310 4.6e-05 ***
uselimit -0.09514 0.10325 -0.921 0.360
useful 0.15932 0.09709 1.641 0.105
T40 0.07096 0.10934 0.649 0.518
outcome 0.06938 0.09025 0.769 0.444
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.54 on 80 degrees of freedom
Multiple R-squared: 0.2848, Adjusted R-squared: 0.2401
F-statistic: 6.373 on 5 and 80 DF, p-value: 4.974e-05
> 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,] 6.770545e-95 1.354109e-94 1.0000000
[2,] 3.087286e-64 6.174571e-64 1.0000000
[3,] 8.184223e-79 1.636845e-78 1.0000000
[4,] 8.505307e-95 1.701061e-94 1.0000000
[5,] 2.759896e-02 5.519793e-02 0.9724010
[6,] 1.213132e-02 2.426265e-02 0.9878687
[7,] 2.033087e-02 4.066174e-02 0.9796691
[8,] 1.096723e-02 2.193445e-02 0.9890328
[9,] 4.961827e-03 9.923654e-03 0.9950382
[10,] 2.245372e-03 4.490744e-03 0.9977546
[11,] 9.785574e-04 1.957115e-03 0.9990214
[12,] 4.839835e-04 9.679671e-04 0.9995160
[13,] 1.729866e-03 3.459731e-03 0.9982701
[14,] 9.432916e-03 1.886583e-02 0.9905671
[15,] 4.057158e-02 8.114315e-02 0.9594284
[16,] 4.102777e-02 8.205554e-02 0.9589722
[17,] 1.040791e-01 2.081582e-01 0.8959209
[18,] 1.212570e-01 2.425141e-01 0.8787430
[19,] 1.024183e-01 2.048366e-01 0.8975817
[20,] 3.083598e-01 6.167196e-01 0.6916402
[21,] 2.578526e-01 5.157053e-01 0.7421474
[22,] 3.107186e-01 6.214372e-01 0.6892814
[23,] 2.567427e-01 5.134853e-01 0.7432573
[24,] 2.151699e-01 4.303397e-01 0.7848301
[25,] 1.880246e-01 3.760492e-01 0.8119754
[26,] 1.786220e-01 3.572440e-01 0.8213780
[27,] 1.418870e-01 2.837740e-01 0.8581130
[28,] 1.109561e-01 2.219121e-01 0.8890439
[29,] 1.447869e-01 2.895738e-01 0.8552131
[30,] 2.976963e-01 5.953927e-01 0.7023037
[31,] 3.230770e-01 6.461541e-01 0.6769230
[32,] 3.629097e-01 7.258194e-01 0.6370903
[33,] 3.025598e-01 6.051197e-01 0.6974402
[34,] 4.455997e-01 8.911994e-01 0.5544003
[35,] 4.198912e-01 8.397824e-01 0.5801088
[36,] 3.704911e-01 7.409821e-01 0.6295089
[37,] 3.524159e-01 7.048317e-01 0.6475841
[38,] 3.529137e-01 7.058273e-01 0.6470863
[39,] 3.037122e-01 6.074245e-01 0.6962878
[40,] 2.640014e-01 5.280028e-01 0.7359986
[41,] 2.654138e-01 5.308276e-01 0.7345862
[42,] 2.234485e-01 4.468970e-01 0.7765515
[43,] 3.886530e-01 7.773060e-01 0.6113470
[44,] 3.306034e-01 6.612068e-01 0.6693966
[45,] 2.912762e-01 5.825525e-01 0.7087238
[46,] 2.444189e-01 4.888377e-01 0.7555811
[47,] 2.008338e-01 4.016676e-01 0.7991662
[48,] 4.670962e-01 9.341925e-01 0.5329038
[49,] 4.878003e-01 9.756006e-01 0.5121997
[50,] 4.370038e-01 8.740077e-01 0.5629962
[51,] 3.891702e-01 7.783404e-01 0.6108298
[52,] 3.591429e-01 7.182857e-01 0.6408571
[53,] 2.979627e-01 5.959254e-01 0.7020373
[54,] 3.614474e-01 7.228948e-01 0.6385526
[55,] 2.967331e-01 5.934662e-01 0.7032669
[56,] 2.462723e-01 4.925446e-01 0.7537277
[57,] 1.917343e-01 3.834686e-01 0.8082657
[58,] 1.452030e-01 2.904060e-01 0.8547970
[59,] 1.057269e-01 2.114539e-01 0.8942731
[60,] 1.134291e-01 2.268583e-01 0.8865709
[61,] 9.021946e-02 1.804389e-01 0.9097805
[62,] 2.993644e-01 5.987289e-01 0.7006356
[63,] 2.218990e-01 4.437980e-01 0.7781010
[64,] 1.805618e-01 3.611236e-01 0.8194382
[65,] 4.047873e-01 8.095747e-01 0.5952127
[66,] 5.204308e-01 9.591384e-01 0.4795692
[67,] 4.110507e-01 8.221014e-01 0.5889493
[68,] 3.221591e-01 6.443181e-01 0.6778409
[69,] 2.531740e-01 5.063481e-01 0.7468260
> postscript(file="/var/fisher/rcomp/tmp/1b98c1355677534.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/fisher/rcomp/tmp/2darg1355677534.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/fisher/rcomp/tmp/3ny5c1355677534.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/fisher/rcomp/tmp/40y8x1355677534.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/fisher/rcomp/tmp/58o4a1355677534.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 = 86
Frequency = 1
1 2 3 4 5 6
-11.0068000 -8.7468917 -8.7468917 -8.7468917 -8.7468917 -15.4250903
7 8 9 10 11 12
-8.7468917 -12.2948012 -12.2158323 -3.9899499 -7.5378594 -8.7468917
13 14 15 16 17 18
33.2869085 -7.5378594 29.8179679 26.2700584 0.8750095 -7.5378594
19 20 21 22 23 24
-12.2158323 -7.3508729 -11.9561497 34.5749097 -20.1820321 -15.4250903
25 26 27 28 29 30
34.2362582 33.2869085 -7.4588905 41.2531083 -12.2158323 -16.7130915
31 32 33 34 35 36
-8.7468917 -3.9899499 -11.9561497 -15.7637418 -8.7468917 -8.7468917
37 38 39 40 41 42
34.4959408 37.7841677 -20.1820321 -20.2610010 -3.8029634 37.7841677
43 44 45 46 47 48
-15.4250903 -7.5378594 -16.7130915 -20.1820321 -8.7468917 -12.2158323
49 50 51 52 53 54
-20.1820321 -8.7468917 37.7051988 0.8750095 -12.2158323 7.6321770
55 56 57 58 59 60
-8.7468917 34.2362582 29.8179679 -12.2158323 -12.2158323 -2.5939311
61 62 63 64 65 66
-11.0068000 33.2869085 -8.7468917 -11.0068000 -8.7468917 -8.7468917
67 68 69 70 71 72
-3.8819323 -3.9899499 -12.2158323 41.2531083 -8.7468917 -12.2158323
73 74 75 76 77 78
37.7841677 46.0100501 -12.2158323 -23.7299416 -12.2158323 29.8179679
79 80 81 82 83 84
0.6153269 -20.2610010 -8.7468917 42.5411095 -8.7468917 7.6321770
85 86
-20.1820321 -3.9899499
> postscript(file="/var/fisher/rcomp/tmp/6kxh01355677534.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 = 86
Frequency = 1
lag(myerror, k = 1) myerror
0 -11.0068000 NA
1 -8.7468917 -11.0068000
2 -8.7468917 -8.7468917
3 -8.7468917 -8.7468917
4 -8.7468917 -8.7468917
5 -15.4250903 -8.7468917
6 -8.7468917 -15.4250903
7 -12.2948012 -8.7468917
8 -12.2158323 -12.2948012
9 -3.9899499 -12.2158323
10 -7.5378594 -3.9899499
11 -8.7468917 -7.5378594
12 33.2869085 -8.7468917
13 -7.5378594 33.2869085
14 29.8179679 -7.5378594
15 26.2700584 29.8179679
16 0.8750095 26.2700584
17 -7.5378594 0.8750095
18 -12.2158323 -7.5378594
19 -7.3508729 -12.2158323
20 -11.9561497 -7.3508729
21 34.5749097 -11.9561497
22 -20.1820321 34.5749097
23 -15.4250903 -20.1820321
24 34.2362582 -15.4250903
25 33.2869085 34.2362582
26 -7.4588905 33.2869085
27 41.2531083 -7.4588905
28 -12.2158323 41.2531083
29 -16.7130915 -12.2158323
30 -8.7468917 -16.7130915
31 -3.9899499 -8.7468917
32 -11.9561497 -3.9899499
33 -15.7637418 -11.9561497
34 -8.7468917 -15.7637418
35 -8.7468917 -8.7468917
36 34.4959408 -8.7468917
37 37.7841677 34.4959408
38 -20.1820321 37.7841677
39 -20.2610010 -20.1820321
40 -3.8029634 -20.2610010
41 37.7841677 -3.8029634
42 -15.4250903 37.7841677
43 -7.5378594 -15.4250903
44 -16.7130915 -7.5378594
45 -20.1820321 -16.7130915
46 -8.7468917 -20.1820321
47 -12.2158323 -8.7468917
48 -20.1820321 -12.2158323
49 -8.7468917 -20.1820321
50 37.7051988 -8.7468917
51 0.8750095 37.7051988
52 -12.2158323 0.8750095
53 7.6321770 -12.2158323
54 -8.7468917 7.6321770
55 34.2362582 -8.7468917
56 29.8179679 34.2362582
57 -12.2158323 29.8179679
58 -12.2158323 -12.2158323
59 -2.5939311 -12.2158323
60 -11.0068000 -2.5939311
61 33.2869085 -11.0068000
62 -8.7468917 33.2869085
63 -11.0068000 -8.7468917
64 -8.7468917 -11.0068000
65 -8.7468917 -8.7468917
66 -3.8819323 -8.7468917
67 -3.9899499 -3.8819323
68 -12.2158323 -3.9899499
69 41.2531083 -12.2158323
70 -8.7468917 41.2531083
71 -12.2158323 -8.7468917
72 37.7841677 -12.2158323
73 46.0100501 37.7841677
74 -12.2158323 46.0100501
75 -23.7299416 -12.2158323
76 -12.2158323 -23.7299416
77 29.8179679 -12.2158323
78 0.6153269 29.8179679
79 -20.2610010 0.6153269
80 -8.7468917 -20.2610010
81 42.5411095 -8.7468917
82 -8.7468917 42.5411095
83 7.6321770 -8.7468917
84 -20.1820321 7.6321770
85 -3.9899499 -20.1820321
86 NA -3.9899499
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -8.7468917 -11.0068000
[2,] -8.7468917 -8.7468917
[3,] -8.7468917 -8.7468917
[4,] -8.7468917 -8.7468917
[5,] -15.4250903 -8.7468917
[6,] -8.7468917 -15.4250903
[7,] -12.2948012 -8.7468917
[8,] -12.2158323 -12.2948012
[9,] -3.9899499 -12.2158323
[10,] -7.5378594 -3.9899499
[11,] -8.7468917 -7.5378594
[12,] 33.2869085 -8.7468917
[13,] -7.5378594 33.2869085
[14,] 29.8179679 -7.5378594
[15,] 26.2700584 29.8179679
[16,] 0.8750095 26.2700584
[17,] -7.5378594 0.8750095
[18,] -12.2158323 -7.5378594
[19,] -7.3508729 -12.2158323
[20,] -11.9561497 -7.3508729
[21,] 34.5749097 -11.9561497
[22,] -20.1820321 34.5749097
[23,] -15.4250903 -20.1820321
[24,] 34.2362582 -15.4250903
[25,] 33.2869085 34.2362582
[26,] -7.4588905 33.2869085
[27,] 41.2531083 -7.4588905
[28,] -12.2158323 41.2531083
[29,] -16.7130915 -12.2158323
[30,] -8.7468917 -16.7130915
[31,] -3.9899499 -8.7468917
[32,] -11.9561497 -3.9899499
[33,] -15.7637418 -11.9561497
[34,] -8.7468917 -15.7637418
[35,] -8.7468917 -8.7468917
[36,] 34.4959408 -8.7468917
[37,] 37.7841677 34.4959408
[38,] -20.1820321 37.7841677
[39,] -20.2610010 -20.1820321
[40,] -3.8029634 -20.2610010
[41,] 37.7841677 -3.8029634
[42,] -15.4250903 37.7841677
[43,] -7.5378594 -15.4250903
[44,] -16.7130915 -7.5378594
[45,] -20.1820321 -16.7130915
[46,] -8.7468917 -20.1820321
[47,] -12.2158323 -8.7468917
[48,] -20.1820321 -12.2158323
[49,] -8.7468917 -20.1820321
[50,] 37.7051988 -8.7468917
[51,] 0.8750095 37.7051988
[52,] -12.2158323 0.8750095
[53,] 7.6321770 -12.2158323
[54,] -8.7468917 7.6321770
[55,] 34.2362582 -8.7468917
[56,] 29.8179679 34.2362582
[57,] -12.2158323 29.8179679
[58,] -12.2158323 -12.2158323
[59,] -2.5939311 -12.2158323
[60,] -11.0068000 -2.5939311
[61,] 33.2869085 -11.0068000
[62,] -8.7468917 33.2869085
[63,] -11.0068000 -8.7468917
[64,] -8.7468917 -11.0068000
[65,] -8.7468917 -8.7468917
[66,] -3.8819323 -8.7468917
[67,] -3.9899499 -3.8819323
[68,] -12.2158323 -3.9899499
[69,] 41.2531083 -12.2158323
[70,] -8.7468917 41.2531083
[71,] -12.2158323 -8.7468917
[72,] 37.7841677 -12.2158323
[73,] 46.0100501 37.7841677
[74,] -12.2158323 46.0100501
[75,] -23.7299416 -12.2158323
[76,] -12.2158323 -23.7299416
[77,] 29.8179679 -12.2158323
[78,] 0.6153269 29.8179679
[79,] -20.2610010 0.6153269
[80,] -8.7468917 -20.2610010
[81,] 42.5411095 -8.7468917
[82,] -8.7468917 42.5411095
[83,] 7.6321770 -8.7468917
[84,] -20.1820321 7.6321770
[85,] -3.9899499 -20.1820321
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -8.7468917 -11.0068000
2 -8.7468917 -8.7468917
3 -8.7468917 -8.7468917
4 -8.7468917 -8.7468917
5 -15.4250903 -8.7468917
6 -8.7468917 -15.4250903
7 -12.2948012 -8.7468917
8 -12.2158323 -12.2948012
9 -3.9899499 -12.2158323
10 -7.5378594 -3.9899499
11 -8.7468917 -7.5378594
12 33.2869085 -8.7468917
13 -7.5378594 33.2869085
14 29.8179679 -7.5378594
15 26.2700584 29.8179679
16 0.8750095 26.2700584
17 -7.5378594 0.8750095
18 -12.2158323 -7.5378594
19 -7.3508729 -12.2158323
20 -11.9561497 -7.3508729
21 34.5749097 -11.9561497
22 -20.1820321 34.5749097
23 -15.4250903 -20.1820321
24 34.2362582 -15.4250903
25 33.2869085 34.2362582
26 -7.4588905 33.2869085
27 41.2531083 -7.4588905
28 -12.2158323 41.2531083
29 -16.7130915 -12.2158323
30 -8.7468917 -16.7130915
31 -3.9899499 -8.7468917
32 -11.9561497 -3.9899499
33 -15.7637418 -11.9561497
34 -8.7468917 -15.7637418
35 -8.7468917 -8.7468917
36 34.4959408 -8.7468917
37 37.7841677 34.4959408
38 -20.1820321 37.7841677
39 -20.2610010 -20.1820321
40 -3.8029634 -20.2610010
41 37.7841677 -3.8029634
42 -15.4250903 37.7841677
43 -7.5378594 -15.4250903
44 -16.7130915 -7.5378594
45 -20.1820321 -16.7130915
46 -8.7468917 -20.1820321
47 -12.2158323 -8.7468917
48 -20.1820321 -12.2158323
49 -8.7468917 -20.1820321
50 37.7051988 -8.7468917
51 0.8750095 37.7051988
52 -12.2158323 0.8750095
53 7.6321770 -12.2158323
54 -8.7468917 7.6321770
55 34.2362582 -8.7468917
56 29.8179679 34.2362582
57 -12.2158323 29.8179679
58 -12.2158323 -12.2158323
59 -2.5939311 -12.2158323
60 -11.0068000 -2.5939311
61 33.2869085 -11.0068000
62 -8.7468917 33.2869085
63 -11.0068000 -8.7468917
64 -8.7468917 -11.0068000
65 -8.7468917 -8.7468917
66 -3.8819323 -8.7468917
67 -3.9899499 -3.8819323
68 -12.2158323 -3.9899499
69 41.2531083 -12.2158323
70 -8.7468917 41.2531083
71 -12.2158323 -8.7468917
72 37.7841677 -12.2158323
73 46.0100501 37.7841677
74 -12.2158323 46.0100501
75 -23.7299416 -12.2158323
76 -12.2158323 -23.7299416
77 29.8179679 -12.2158323
78 0.6153269 29.8179679
79 -20.2610010 0.6153269
80 -8.7468917 -20.2610010
81 42.5411095 -8.7468917
82 -8.7468917 42.5411095
83 7.6321770 -8.7468917
84 -20.1820321 7.6321770
85 -3.9899499 -20.1820321
> 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/fisher/rcomp/tmp/7ui1d1355677534.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/fisher/rcomp/tmp/8vzlr1355677534.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/fisher/rcomp/tmp/972oz1355677534.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/fisher/rcomp/tmp/10se7a1355677534.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/fisher/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/fisher/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/fisher/rcomp/tmp/11ich61355677534.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/fisher/rcomp/tmp/12b4nr1355677534.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/fisher/rcomp/tmp/13f9ic1355677534.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/fisher/rcomp/tmp/14i8eh1355677534.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/fisher/rcomp/tmp/153kix1355677534.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/fisher/rcomp/tmp/16vunu1355677534.tab")
+ }
>
> try(system("convert tmp/1b98c1355677534.ps tmp/1b98c1355677534.png",intern=TRUE))
character(0)
> try(system("convert tmp/2darg1355677534.ps tmp/2darg1355677534.png",intern=TRUE))
character(0)
> try(system("convert tmp/3ny5c1355677534.ps tmp/3ny5c1355677534.png",intern=TRUE))
character(0)
> try(system("convert tmp/40y8x1355677534.ps tmp/40y8x1355677534.png",intern=TRUE))
character(0)
> try(system("convert tmp/58o4a1355677534.ps tmp/58o4a1355677534.png",intern=TRUE))
character(0)
> try(system("convert tmp/6kxh01355677534.ps tmp/6kxh01355677534.png",intern=TRUE))
character(0)
> try(system("convert tmp/7ui1d1355677534.ps tmp/7ui1d1355677534.png",intern=TRUE))
character(0)
> try(system("convert tmp/8vzlr1355677534.ps tmp/8vzlr1355677534.png",intern=TRUE))
character(0)
> try(system("convert tmp/972oz1355677534.ps tmp/972oz1355677534.png",intern=TRUE))
character(0)
> try(system("convert tmp/10se7a1355677534.ps tmp/10se7a1355677534.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.853 1.717 8.569