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.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(101.02,0,100.67,0,100.47,0,100.38,0,100.33,0,100.34,0,100.37,0,100.39,0,100.21,0,100.21,0,100.22,0,100.28,0,100.25,0,100.25,0,100.21,0,100.16,0,100.18,0,100.1,1,99.96,1,99.88,1,99.88,1,99.86,1,99.84,1,99.8,1,99.82,1,99.81,1,99.92,1,100.03,1,99.99,1,100.02,1,100.01,1,100.13,1,100.33,1,100.13,1,99.96,1,100.05,1,99.83,1,99.8,1,100.01,1,100.1,1,100.13,1,100.16,1,100.41,1,101.34,1,101.65,1,101.85,1,102.07,1,102.12,1,102.14,1,102.21,1,102.28,1,102.19,1,102.33,1,102.54,1,102.44,1,102.78,1,102.9,1,103.08,1,102.77,1,102.65,1,102.71,1,103.29,1,102.86,1,103.45,1,103.72,1,103.65,1,103.83,1,104.45,1,105.14,1,105.07,1,105.31,1,105.19,1,105.3,1,105.02,1,105.17,1,105.28,1,105.45,1,105.38,1,105.8,1,105.96,1,105.08,1,105.11,1,105.61,1,105.5,1),dim=c(2,84),dimnames=list(c('Suiker','Dummy'),1:84))
> y <- array(NA,dim=c(2,84),dimnames=list(c('Suiker','Dummy'),1:84))
> 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
Suiker Dummy M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 101.02 0 1 0 0 0 0 0 0 0 0 0 0
2 100.67 0 0 1 0 0 0 0 0 0 0 0 0
3 100.47 0 0 0 1 0 0 0 0 0 0 0 0
4 100.38 0 0 0 0 1 0 0 0 0 0 0 0
5 100.33 0 0 0 0 0 1 0 0 0 0 0 0
6 100.34 0 0 0 0 0 0 1 0 0 0 0 0
7 100.37 0 0 0 0 0 0 0 1 0 0 0 0
8 100.39 0 0 0 0 0 0 0 0 1 0 0 0
9 100.21 0 0 0 0 0 0 0 0 0 1 0 0
10 100.21 0 0 0 0 0 0 0 0 0 0 1 0
11 100.22 0 0 0 0 0 0 0 0 0 0 0 1
12 100.28 0 0 0 0 0 0 0 0 0 0 0 0
13 100.25 0 1 0 0 0 0 0 0 0 0 0 0
14 100.25 0 0 1 0 0 0 0 0 0 0 0 0
15 100.21 0 0 0 1 0 0 0 0 0 0 0 0
16 100.16 0 0 0 0 1 0 0 0 0 0 0 0
17 100.18 0 0 0 0 0 1 0 0 0 0 0 0
18 100.10 1 0 0 0 0 0 1 0 0 0 0 0
19 99.96 1 0 0 0 0 0 0 1 0 0 0 0
20 99.88 1 0 0 0 0 0 0 0 1 0 0 0
21 99.88 1 0 0 0 0 0 0 0 0 1 0 0
22 99.86 1 0 0 0 0 0 0 0 0 0 1 0
23 99.84 1 0 0 0 0 0 0 0 0 0 0 1
24 99.80 1 0 0 0 0 0 0 0 0 0 0 0
25 99.82 1 1 0 0 0 0 0 0 0 0 0 0
26 99.81 1 0 1 0 0 0 0 0 0 0 0 0
27 99.92 1 0 0 1 0 0 0 0 0 0 0 0
28 100.03 1 0 0 0 1 0 0 0 0 0 0 0
29 99.99 1 0 0 0 0 1 0 0 0 0 0 0
30 100.02 1 0 0 0 0 0 1 0 0 0 0 0
31 100.01 1 0 0 0 0 0 0 1 0 0 0 0
32 100.13 1 0 0 0 0 0 0 0 1 0 0 0
33 100.33 1 0 0 0 0 0 0 0 0 1 0 0
34 100.13 1 0 0 0 0 0 0 0 0 0 1 0
35 99.96 1 0 0 0 0 0 0 0 0 0 0 1
36 100.05 1 0 0 0 0 0 0 0 0 0 0 0
37 99.83 1 1 0 0 0 0 0 0 0 0 0 0
38 99.80 1 0 1 0 0 0 0 0 0 0 0 0
39 100.01 1 0 0 1 0 0 0 0 0 0 0 0
40 100.10 1 0 0 0 1 0 0 0 0 0 0 0
41 100.13 1 0 0 0 0 1 0 0 0 0 0 0
42 100.16 1 0 0 0 0 0 1 0 0 0 0 0
43 100.41 1 0 0 0 0 0 0 1 0 0 0 0
44 101.34 1 0 0 0 0 0 0 0 1 0 0 0
45 101.65 1 0 0 0 0 0 0 0 0 1 0 0
46 101.85 1 0 0 0 0 0 0 0 0 0 1 0
47 102.07 1 0 0 0 0 0 0 0 0 0 0 1
48 102.12 1 0 0 0 0 0 0 0 0 0 0 0
49 102.14 1 1 0 0 0 0 0 0 0 0 0 0
50 102.21 1 0 1 0 0 0 0 0 0 0 0 0
51 102.28 1 0 0 1 0 0 0 0 0 0 0 0
52 102.19 1 0 0 0 1 0 0 0 0 0 0 0
53 102.33 1 0 0 0 0 1 0 0 0 0 0 0
54 102.54 1 0 0 0 0 0 1 0 0 0 0 0
55 102.44 1 0 0 0 0 0 0 1 0 0 0 0
56 102.78 1 0 0 0 0 0 0 0 1 0 0 0
57 102.90 1 0 0 0 0 0 0 0 0 1 0 0
58 103.08 1 0 0 0 0 0 0 0 0 0 1 0
59 102.77 1 0 0 0 0 0 0 0 0 0 0 1
60 102.65 1 0 0 0 0 0 0 0 0 0 0 0
61 102.71 1 1 0 0 0 0 0 0 0 0 0 0
62 103.29 1 0 1 0 0 0 0 0 0 0 0 0
63 102.86 1 0 0 1 0 0 0 0 0 0 0 0
64 103.45 1 0 0 0 1 0 0 0 0 0 0 0
65 103.72 1 0 0 0 0 1 0 0 0 0 0 0
66 103.65 1 0 0 0 0 0 1 0 0 0 0 0
67 103.83 1 0 0 0 0 0 0 1 0 0 0 0
68 104.45 1 0 0 0 0 0 0 0 1 0 0 0
69 105.14 1 0 0 0 0 0 0 0 0 1 0 0
70 105.07 1 0 0 0 0 0 0 0 0 0 1 0
71 105.31 1 0 0 0 0 0 0 0 0 0 0 1
72 105.19 1 0 0 0 0 0 0 0 0 0 0 0
73 105.30 1 1 0 0 0 0 0 0 0 0 0 0
74 105.02 1 0 1 0 0 0 0 0 0 0 0 0
75 105.17 1 0 0 1 0 0 0 0 0 0 0 0
76 105.28 1 0 0 0 1 0 0 0 0 0 0 0
77 105.45 1 0 0 0 0 1 0 0 0 0 0 0
78 105.38 1 0 0 0 0 0 1 0 0 0 0 0
79 105.80 1 0 0 0 0 0 0 1 0 0 0 0
80 105.96 1 0 0 0 0 0 0 0 1 0 0 0
81 105.08 1 0 0 0 0 0 0 0 0 1 0 0
82 105.11 1 0 0 0 0 0 0 0 0 0 1 0
83 105.61 1 0 0 0 0 0 0 0 0 0 0 1
84 105.50 1 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) Dummy M1 M2 M3 M4
100.60835 1.88859 -0.37592 -0.37877 -0.39734 -0.30163
M5 M6 M7 M8 M9 M10
-0.22449 -0.48571 -0.39571 -0.09429 -0.05714 -0.04000
M11
0.02714
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.69694 -2.09000 0.01827 1.19255 3.69877
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 100.60835 0.91932 109.438 < 2e-16 ***
Dummy 1.88859 0.56880 3.320 0.00142 **
M1 -0.37592 1.10522 -0.340 0.73476
M2 -0.37877 1.10522 -0.343 0.73283
M3 -0.39734 1.10522 -0.360 0.72028
M4 -0.30163 1.10522 -0.273 0.78571
M5 -0.22449 1.10522 -0.203 0.83963
M6 -0.48571 1.10223 -0.441 0.66079
M7 -0.39571 1.10223 -0.359 0.72065
M8 -0.09429 1.10223 -0.086 0.93207
M9 -0.05714 1.10223 -0.052 0.95880
M10 -0.04000 1.10223 -0.036 0.97115
M11 0.02714 1.10223 0.025 0.98042
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.062 on 71 degrees of freedom
Multiple R-squared: 0.1494, Adjusted R-squared: 0.005631
F-statistic: 1.039 on 12 and 71 DF, p-value: 0.4239
> 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.073370e-03 1.214674e-02 0.993926630
[2,] 8.402222e-04 1.680444e-03 0.999159778
[3,] 1.041326e-04 2.082651e-04 0.999895867
[4,] 1.266542e-05 2.533083e-05 0.999987335
[5,] 1.561667e-06 3.123334e-06 0.999998438
[6,] 1.745934e-07 3.491868e-07 0.999999825
[7,] 1.885350e-08 3.770701e-08 0.999999981
[8,] 2.009853e-09 4.019705e-09 0.999999998
[9,] 2.239302e-10 4.478604e-10 1.000000000
[10,] 6.878295e-11 1.375659e-10 1.000000000
[11,] 9.164401e-12 1.832880e-11 1.000000000
[12,] 9.741428e-13 1.948286e-12 1.000000000
[13,] 1.481038e-13 2.962076e-13 1.000000000
[14,] 2.048768e-14 4.097536e-14 1.000000000
[15,] 2.243640e-15 4.487279e-15 1.000000000
[16,] 2.640721e-16 5.281443e-16 1.000000000
[17,] 5.008762e-17 1.001752e-16 1.000000000
[18,] 5.055100e-17 1.011020e-16 1.000000000
[19,] 1.606109e-17 3.212218e-17 1.000000000
[20,] 3.721257e-18 7.442514e-18 1.000000000
[21,] 1.092779e-18 2.185559e-18 1.000000000
[22,] 4.533605e-19 9.067211e-19 1.000000000
[23,] 1.500599e-19 3.001199e-19 1.000000000
[24,] 3.766224e-20 7.532447e-20 1.000000000
[25,] 1.411112e-20 2.822225e-20 1.000000000
[26,] 8.272912e-21 1.654582e-20 1.000000000
[27,] 4.545691e-21 9.091383e-21 1.000000000
[28,] 1.316307e-20 2.632614e-20 1.000000000
[29,] 2.759252e-16 5.518505e-16 1.000000000
[30,] 6.390426e-13 1.278085e-12 1.000000000
[31,] 3.012414e-10 6.024828e-10 1.000000000
[32,] 4.849102e-08 9.698204e-08 0.999999952
[33,] 1.406611e-06 2.813223e-06 0.999998593
[34,] 9.650558e-06 1.930112e-05 0.999990349
[35,] 5.629891e-05 1.125978e-04 0.999943701
[36,] 2.049253e-04 4.098506e-04 0.999795075
[37,] 6.167156e-04 1.233431e-03 0.999383284
[38,] 1.859764e-03 3.719527e-03 0.998140236
[39,] 4.995129e-03 9.990258e-03 0.995004871
[40,] 1.269241e-02 2.538482e-02 0.987307591
[41,] 3.070304e-02 6.140608e-02 0.969296961
[42,] 5.919222e-02 1.183844e-01 0.940807783
[43,] 9.975885e-02 1.995177e-01 0.900241146
[44,] 2.003553e-01 4.007106e-01 0.799644710
[45,] 3.646505e-01 7.293011e-01 0.635349464
[46,] 5.063645e-01 9.872711e-01 0.493635544
[47,] 5.661350e-01 8.677300e-01 0.433864984
[48,] 6.888803e-01 6.222395e-01 0.311119748
[49,] 7.565915e-01 4.868170e-01 0.243408484
[50,] 8.101690e-01 3.796620e-01 0.189831014
[51,] 8.607790e-01 2.784420e-01 0.139221011
[52,] 9.512202e-01 9.755965e-02 0.048779826
[53,] 9.970544e-01 5.891201e-03 0.002945601
> postscript(file="/var/www/html/freestat/rcomp/tmp/16vch1229365756.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/freestat/rcomp/tmp/284lg1229365756.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/freestat/rcomp/tmp/3ie1o1229365756.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/freestat/rcomp/tmp/4qybl1229365756.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/freestat/rcomp/tmp/5pkm71229365756.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 = 84
Frequency = 1
1 2 3 4 5 6
0.787562112 0.440419255 0.258990683 0.073276398 -0.053866460 0.217360248
7 8 9 10 11 12
0.157360248 -0.124068323 -0.341211180 -0.358354037 -0.415496894 -0.328354037
13 14 15 16 17 18
0.017562112 0.020419255 -0.001009317 -0.146723602 -0.203866460 -1.911226708
19 20 21 22 23 24
-2.141226708 -2.522655280 -2.559798137 -2.596940994 -2.684083851 -2.696940994
25 26 27 28 29 30
-2.301024845 -2.308167702 -2.179596273 -2.165310559 -2.282453416 -1.991226708
31 32 33 34 35 36
-2.091226708 -2.272655280 -2.109798137 -2.326940994 -2.564083851 -2.446940994
37 38 39 40 41 42
-2.291024845 -2.318167702 -2.089596273 -2.095310559 -2.142453416 -1.851226708
43 44 45 46 47 48
-1.691226708 -1.062655280 -0.789798137 -0.606940994 -0.454083851 -0.376940994
49 50 51 52 53 54
0.018975155 0.091832298 0.180403727 -0.005310559 0.057546584 0.528773292
55 56 57 58 59 60
0.338773292 0.377344720 0.460201863 0.623059006 0.245916149 0.153059006
61 62 63 64 65 66
0.588975155 1.171832298 0.760403727 1.254689441 1.447546584 1.638773292
67 68 69 70 71 72
1.728773292 2.047344720 2.700201863 2.613059006 2.785916149 2.693059006
73 74 75 76 77 78
3.178975155 2.901832298 3.070403727 3.084689441 3.177546584 3.368773292
79 80 81 82 83 84
3.698773292 3.557344720 2.640201863 2.653059006 3.085916149 3.003059006
> postscript(file="/var/www/html/freestat/rcomp/tmp/6wikf1229365756.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 = 84
Frequency = 1
lag(myerror, k = 1) myerror
0 0.787562112 NA
1 0.440419255 0.787562112
2 0.258990683 0.440419255
3 0.073276398 0.258990683
4 -0.053866460 0.073276398
5 0.217360248 -0.053866460
6 0.157360248 0.217360248
7 -0.124068323 0.157360248
8 -0.341211180 -0.124068323
9 -0.358354037 -0.341211180
10 -0.415496894 -0.358354037
11 -0.328354037 -0.415496894
12 0.017562112 -0.328354037
13 0.020419255 0.017562112
14 -0.001009317 0.020419255
15 -0.146723602 -0.001009317
16 -0.203866460 -0.146723602
17 -1.911226708 -0.203866460
18 -2.141226708 -1.911226708
19 -2.522655280 -2.141226708
20 -2.559798137 -2.522655280
21 -2.596940994 -2.559798137
22 -2.684083851 -2.596940994
23 -2.696940994 -2.684083851
24 -2.301024845 -2.696940994
25 -2.308167702 -2.301024845
26 -2.179596273 -2.308167702
27 -2.165310559 -2.179596273
28 -2.282453416 -2.165310559
29 -1.991226708 -2.282453416
30 -2.091226708 -1.991226708
31 -2.272655280 -2.091226708
32 -2.109798137 -2.272655280
33 -2.326940994 -2.109798137
34 -2.564083851 -2.326940994
35 -2.446940994 -2.564083851
36 -2.291024845 -2.446940994
37 -2.318167702 -2.291024845
38 -2.089596273 -2.318167702
39 -2.095310559 -2.089596273
40 -2.142453416 -2.095310559
41 -1.851226708 -2.142453416
42 -1.691226708 -1.851226708
43 -1.062655280 -1.691226708
44 -0.789798137 -1.062655280
45 -0.606940994 -0.789798137
46 -0.454083851 -0.606940994
47 -0.376940994 -0.454083851
48 0.018975155 -0.376940994
49 0.091832298 0.018975155
50 0.180403727 0.091832298
51 -0.005310559 0.180403727
52 0.057546584 -0.005310559
53 0.528773292 0.057546584
54 0.338773292 0.528773292
55 0.377344720 0.338773292
56 0.460201863 0.377344720
57 0.623059006 0.460201863
58 0.245916149 0.623059006
59 0.153059006 0.245916149
60 0.588975155 0.153059006
61 1.171832298 0.588975155
62 0.760403727 1.171832298
63 1.254689441 0.760403727
64 1.447546584 1.254689441
65 1.638773292 1.447546584
66 1.728773292 1.638773292
67 2.047344720 1.728773292
68 2.700201863 2.047344720
69 2.613059006 2.700201863
70 2.785916149 2.613059006
71 2.693059006 2.785916149
72 3.178975155 2.693059006
73 2.901832298 3.178975155
74 3.070403727 2.901832298
75 3.084689441 3.070403727
76 3.177546584 3.084689441
77 3.368773292 3.177546584
78 3.698773292 3.368773292
79 3.557344720 3.698773292
80 2.640201863 3.557344720
81 2.653059006 2.640201863
82 3.085916149 2.653059006
83 3.003059006 3.085916149
84 NA 3.003059006
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.440419255 0.787562112
[2,] 0.258990683 0.440419255
[3,] 0.073276398 0.258990683
[4,] -0.053866460 0.073276398
[5,] 0.217360248 -0.053866460
[6,] 0.157360248 0.217360248
[7,] -0.124068323 0.157360248
[8,] -0.341211180 -0.124068323
[9,] -0.358354037 -0.341211180
[10,] -0.415496894 -0.358354037
[11,] -0.328354037 -0.415496894
[12,] 0.017562112 -0.328354037
[13,] 0.020419255 0.017562112
[14,] -0.001009317 0.020419255
[15,] -0.146723602 -0.001009317
[16,] -0.203866460 -0.146723602
[17,] -1.911226708 -0.203866460
[18,] -2.141226708 -1.911226708
[19,] -2.522655280 -2.141226708
[20,] -2.559798137 -2.522655280
[21,] -2.596940994 -2.559798137
[22,] -2.684083851 -2.596940994
[23,] -2.696940994 -2.684083851
[24,] -2.301024845 -2.696940994
[25,] -2.308167702 -2.301024845
[26,] -2.179596273 -2.308167702
[27,] -2.165310559 -2.179596273
[28,] -2.282453416 -2.165310559
[29,] -1.991226708 -2.282453416
[30,] -2.091226708 -1.991226708
[31,] -2.272655280 -2.091226708
[32,] -2.109798137 -2.272655280
[33,] -2.326940994 -2.109798137
[34,] -2.564083851 -2.326940994
[35,] -2.446940994 -2.564083851
[36,] -2.291024845 -2.446940994
[37,] -2.318167702 -2.291024845
[38,] -2.089596273 -2.318167702
[39,] -2.095310559 -2.089596273
[40,] -2.142453416 -2.095310559
[41,] -1.851226708 -2.142453416
[42,] -1.691226708 -1.851226708
[43,] -1.062655280 -1.691226708
[44,] -0.789798137 -1.062655280
[45,] -0.606940994 -0.789798137
[46,] -0.454083851 -0.606940994
[47,] -0.376940994 -0.454083851
[48,] 0.018975155 -0.376940994
[49,] 0.091832298 0.018975155
[50,] 0.180403727 0.091832298
[51,] -0.005310559 0.180403727
[52,] 0.057546584 -0.005310559
[53,] 0.528773292 0.057546584
[54,] 0.338773292 0.528773292
[55,] 0.377344720 0.338773292
[56,] 0.460201863 0.377344720
[57,] 0.623059006 0.460201863
[58,] 0.245916149 0.623059006
[59,] 0.153059006 0.245916149
[60,] 0.588975155 0.153059006
[61,] 1.171832298 0.588975155
[62,] 0.760403727 1.171832298
[63,] 1.254689441 0.760403727
[64,] 1.447546584 1.254689441
[65,] 1.638773292 1.447546584
[66,] 1.728773292 1.638773292
[67,] 2.047344720 1.728773292
[68,] 2.700201863 2.047344720
[69,] 2.613059006 2.700201863
[70,] 2.785916149 2.613059006
[71,] 2.693059006 2.785916149
[72,] 3.178975155 2.693059006
[73,] 2.901832298 3.178975155
[74,] 3.070403727 2.901832298
[75,] 3.084689441 3.070403727
[76,] 3.177546584 3.084689441
[77,] 3.368773292 3.177546584
[78,] 3.698773292 3.368773292
[79,] 3.557344720 3.698773292
[80,] 2.640201863 3.557344720
[81,] 2.653059006 2.640201863
[82,] 3.085916149 2.653059006
[83,] 3.003059006 3.085916149
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.440419255 0.787562112
2 0.258990683 0.440419255
3 0.073276398 0.258990683
4 -0.053866460 0.073276398
5 0.217360248 -0.053866460
6 0.157360248 0.217360248
7 -0.124068323 0.157360248
8 -0.341211180 -0.124068323
9 -0.358354037 -0.341211180
10 -0.415496894 -0.358354037
11 -0.328354037 -0.415496894
12 0.017562112 -0.328354037
13 0.020419255 0.017562112
14 -0.001009317 0.020419255
15 -0.146723602 -0.001009317
16 -0.203866460 -0.146723602
17 -1.911226708 -0.203866460
18 -2.141226708 -1.911226708
19 -2.522655280 -2.141226708
20 -2.559798137 -2.522655280
21 -2.596940994 -2.559798137
22 -2.684083851 -2.596940994
23 -2.696940994 -2.684083851
24 -2.301024845 -2.696940994
25 -2.308167702 -2.301024845
26 -2.179596273 -2.308167702
27 -2.165310559 -2.179596273
28 -2.282453416 -2.165310559
29 -1.991226708 -2.282453416
30 -2.091226708 -1.991226708
31 -2.272655280 -2.091226708
32 -2.109798137 -2.272655280
33 -2.326940994 -2.109798137
34 -2.564083851 -2.326940994
35 -2.446940994 -2.564083851
36 -2.291024845 -2.446940994
37 -2.318167702 -2.291024845
38 -2.089596273 -2.318167702
39 -2.095310559 -2.089596273
40 -2.142453416 -2.095310559
41 -1.851226708 -2.142453416
42 -1.691226708 -1.851226708
43 -1.062655280 -1.691226708
44 -0.789798137 -1.062655280
45 -0.606940994 -0.789798137
46 -0.454083851 -0.606940994
47 -0.376940994 -0.454083851
48 0.018975155 -0.376940994
49 0.091832298 0.018975155
50 0.180403727 0.091832298
51 -0.005310559 0.180403727
52 0.057546584 -0.005310559
53 0.528773292 0.057546584
54 0.338773292 0.528773292
55 0.377344720 0.338773292
56 0.460201863 0.377344720
57 0.623059006 0.460201863
58 0.245916149 0.623059006
59 0.153059006 0.245916149
60 0.588975155 0.153059006
61 1.171832298 0.588975155
62 0.760403727 1.171832298
63 1.254689441 0.760403727
64 1.447546584 1.254689441
65 1.638773292 1.447546584
66 1.728773292 1.638773292
67 2.047344720 1.728773292
68 2.700201863 2.047344720
69 2.613059006 2.700201863
70 2.785916149 2.613059006
71 2.693059006 2.785916149
72 3.178975155 2.693059006
73 2.901832298 3.178975155
74 3.070403727 2.901832298
75 3.084689441 3.070403727
76 3.177546584 3.084689441
77 3.368773292 3.177546584
78 3.698773292 3.368773292
79 3.557344720 3.698773292
80 2.640201863 3.557344720
81 2.653059006 2.640201863
82 3.085916149 2.653059006
83 3.003059006 3.085916149
> 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/7dgoi1229365756.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/freestat/rcomp/tmp/8qi9l1229365756.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/freestat/rcomp/tmp/95ufw1229365756.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/freestat/rcomp/tmp/103lhq1229365756.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/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/11ccbh1229365756.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/12ylot1229365756.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/13ui4u1229365756.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/14m8ph1229365756.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/15c2m81229365756.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/16ke141229365756.tab")
+ }
>
> system("convert tmp/16vch1229365756.ps tmp/16vch1229365756.png")
> system("convert tmp/284lg1229365756.ps tmp/284lg1229365756.png")
> system("convert tmp/3ie1o1229365756.ps tmp/3ie1o1229365756.png")
> system("convert tmp/4qybl1229365756.ps tmp/4qybl1229365756.png")
> system("convert tmp/5pkm71229365756.ps tmp/5pkm71229365756.png")
> system("convert tmp/6wikf1229365756.ps tmp/6wikf1229365756.png")
> system("convert tmp/7dgoi1229365756.ps tmp/7dgoi1229365756.png")
> system("convert tmp/8qi9l1229365756.ps tmp/8qi9l1229365756.png")
> system("convert tmp/95ufw1229365756.ps tmp/95ufw1229365756.png")
> system("convert tmp/103lhq1229365756.ps tmp/103lhq1229365756.png")
>
>
> proc.time()
user system elapsed
4.026 2.520 4.816