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(1.43,0,1.43,0,1.43,0,1.43,0,1.43,0,1.43,0,1.43,0,1.43,0,1.43,0,1.43,0,1.43,0,1.43,0,1.43,0,1.43,0,1.43,0,1.43,0,1.43,0,1.43,0,1.44,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.48,0,1.57,1,1.58,1,1.58,1,1.58,1,1.58,1,1.59,1,1.6,1,1.6,1,1.61,1,1.61,1,1.61,1,1.62,1,1.63,1,1.63,1,1.64,1,1.64,1,1.64,1,1.64,1,1.64,1,1.65,1,1.65,1,1.65,1,1.65,1,1.65,1,1.66,1,1.66,1,1.67,1,1.68,1,1.68,1,1.68,1,1.68,1,1.69,1,1.7,1,1.7,1,1.71,1,1.72,1,1.73,1,1.74,1,1.74,1,1.75,1,1.75,1,1.75,1,1.76,1,1.79,1,1.83,1,1.84,1,1.85,1),dim=c(2,96),dimnames=list(c('Y','X'),1:96))
> y <- array(NA,dim=c(2,96),dimnames=list(c('Y','X'),1:96))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 1.43 0 1 0 0 0 0 0 0 0 0 0 0 1
2 1.43 0 0 1 0 0 0 0 0 0 0 0 0 2
3 1.43 0 0 0 1 0 0 0 0 0 0 0 0 3
4 1.43 0 0 0 0 1 0 0 0 0 0 0 0 4
5 1.43 0 0 0 0 0 1 0 0 0 0 0 0 5
6 1.43 0 0 0 0 0 0 1 0 0 0 0 0 6
7 1.43 0 0 0 0 0 0 0 1 0 0 0 0 7
8 1.43 0 0 0 0 0 0 0 0 1 0 0 0 8
9 1.43 0 0 0 0 0 0 0 0 0 1 0 0 9
10 1.43 0 0 0 0 0 0 0 0 0 0 1 0 10
11 1.43 0 0 0 0 0 0 0 0 0 0 0 1 11
12 1.43 0 0 0 0 0 0 0 0 0 0 0 0 12
13 1.43 0 1 0 0 0 0 0 0 0 0 0 0 13
14 1.43 0 0 1 0 0 0 0 0 0 0 0 0 14
15 1.43 0 0 0 1 0 0 0 0 0 0 0 0 15
16 1.43 0 0 0 0 1 0 0 0 0 0 0 0 16
17 1.43 0 0 0 0 0 1 0 0 0 0 0 0 17
18 1.43 0 0 0 0 0 0 1 0 0 0 0 0 18
19 1.44 0 0 0 0 0 0 0 1 0 0 0 0 19
20 1.48 0 0 0 0 0 0 0 0 1 0 0 0 20
21 1.48 0 0 0 0 0 0 0 0 0 1 0 0 21
22 1.48 0 0 0 0 0 0 0 0 0 0 1 0 22
23 1.48 0 0 0 0 0 0 0 0 0 0 0 1 23
24 1.48 0 0 0 0 0 0 0 0 0 0 0 0 24
25 1.48 0 1 0 0 0 0 0 0 0 0 0 0 25
26 1.48 0 0 1 0 0 0 0 0 0 0 0 0 26
27 1.48 0 0 0 1 0 0 0 0 0 0 0 0 27
28 1.48 0 0 0 0 1 0 0 0 0 0 0 0 28
29 1.48 0 0 0 0 0 1 0 0 0 0 0 0 29
30 1.48 0 0 0 0 0 0 1 0 0 0 0 0 30
31 1.48 0 0 0 0 0 0 0 1 0 0 0 0 31
32 1.48 0 0 0 0 0 0 0 0 1 0 0 0 32
33 1.48 0 0 0 0 0 0 0 0 0 1 0 0 33
34 1.48 0 0 0 0 0 0 0 0 0 0 1 0 34
35 1.48 0 0 0 0 0 0 0 0 0 0 0 1 35
36 1.48 0 0 0 0 0 0 0 0 0 0 0 0 36
37 1.48 0 1 0 0 0 0 0 0 0 0 0 0 37
38 1.48 0 0 1 0 0 0 0 0 0 0 0 0 38
39 1.48 0 0 0 1 0 0 0 0 0 0 0 0 39
40 1.48 0 0 0 0 1 0 0 0 0 0 0 0 40
41 1.48 0 0 0 0 0 1 0 0 0 0 0 0 41
42 1.48 0 0 0 0 0 0 1 0 0 0 0 0 42
43 1.48 0 0 0 0 0 0 0 1 0 0 0 0 43
44 1.48 0 0 0 0 0 0 0 0 1 0 0 0 44
45 1.48 0 0 0 0 0 0 0 0 0 1 0 0 45
46 1.48 0 0 0 0 0 0 0 0 0 0 1 0 46
47 1.48 0 0 0 0 0 0 0 0 0 0 0 1 47
48 1.48 0 0 0 0 0 0 0 0 0 0 0 0 48
49 1.48 0 1 0 0 0 0 0 0 0 0 0 0 49
50 1.57 1 0 1 0 0 0 0 0 0 0 0 0 50
51 1.58 1 0 0 1 0 0 0 0 0 0 0 0 51
52 1.58 1 0 0 0 1 0 0 0 0 0 0 0 52
53 1.58 1 0 0 0 0 1 0 0 0 0 0 0 53
54 1.58 1 0 0 0 0 0 1 0 0 0 0 0 54
55 1.59 1 0 0 0 0 0 0 1 0 0 0 0 55
56 1.60 1 0 0 0 0 0 0 0 1 0 0 0 56
57 1.60 1 0 0 0 0 0 0 0 0 1 0 0 57
58 1.61 1 0 0 0 0 0 0 0 0 0 1 0 58
59 1.61 1 0 0 0 0 0 0 0 0 0 0 1 59
60 1.61 1 0 0 0 0 0 0 0 0 0 0 0 60
61 1.62 1 1 0 0 0 0 0 0 0 0 0 0 61
62 1.63 1 0 1 0 0 0 0 0 0 0 0 0 62
63 1.63 1 0 0 1 0 0 0 0 0 0 0 0 63
64 1.64 1 0 0 0 1 0 0 0 0 0 0 0 64
65 1.64 1 0 0 0 0 1 0 0 0 0 0 0 65
66 1.64 1 0 0 0 0 0 1 0 0 0 0 0 66
67 1.64 1 0 0 0 0 0 0 1 0 0 0 0 67
68 1.64 1 0 0 0 0 0 0 0 1 0 0 0 68
69 1.65 1 0 0 0 0 0 0 0 0 1 0 0 69
70 1.65 1 0 0 0 0 0 0 0 0 0 1 0 70
71 1.65 1 0 0 0 0 0 0 0 0 0 0 1 71
72 1.65 1 0 0 0 0 0 0 0 0 0 0 0 72
73 1.65 1 1 0 0 0 0 0 0 0 0 0 0 73
74 1.66 1 0 1 0 0 0 0 0 0 0 0 0 74
75 1.66 1 0 0 1 0 0 0 0 0 0 0 0 75
76 1.67 1 0 0 0 1 0 0 0 0 0 0 0 76
77 1.68 1 0 0 0 0 1 0 0 0 0 0 0 77
78 1.68 1 0 0 0 0 0 1 0 0 0 0 0 78
79 1.68 1 0 0 0 0 0 0 1 0 0 0 0 79
80 1.68 1 0 0 0 0 0 0 0 1 0 0 0 80
81 1.69 1 0 0 0 0 0 0 0 0 1 0 0 81
82 1.70 1 0 0 0 0 0 0 0 0 0 1 0 82
83 1.70 1 0 0 0 0 0 0 0 0 0 0 1 83
84 1.71 1 0 0 0 0 0 0 0 0 0 0 0 84
85 1.72 1 1 0 0 0 0 0 0 0 0 0 0 85
86 1.73 1 0 1 0 0 0 0 0 0 0 0 0 86
87 1.74 1 0 0 1 0 0 0 0 0 0 0 0 87
88 1.74 1 0 0 0 1 0 0 0 0 0 0 0 88
89 1.75 1 0 0 0 0 1 0 0 0 0 0 0 89
90 1.75 1 0 0 0 0 0 1 0 0 0 0 0 90
91 1.75 1 0 0 0 0 0 0 1 0 0 0 0 91
92 1.76 1 0 0 0 0 0 0 0 1 0 0 0 92
93 1.79 1 0 0 0 0 0 0 0 0 1 0 0 93
94 1.83 1 0 0 0 0 0 0 0 0 0 1 0 94
95 1.84 1 0 0 0 0 0 0 0 0 0 0 1 95
96 1.85 1 0 0 0 0 0 0 0 0 0 0 0 96
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
1.3888686 0.0646399 -0.0082964 -0.0044331 -0.0049898 -0.0055464
M5 M6 M7 M8 M9 M10
-0.0061031 -0.0091598 -0.0097165 -0.0052732 -0.0020799 0.0023634
M11 t
0.0005567 0.0030567
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.05559 -0.02003 -0.00500 0.01784 0.10305
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.3888686 0.0135969 102.146 < 2e-16 ***
X 0.0646399 0.0130898 4.938 4.09e-06 ***
M1 -0.0082964 0.0157992 -0.525 0.601
M2 -0.0044331 0.0159121 -0.279 0.781
M3 -0.0049898 0.0158784 -0.314 0.754
M4 -0.0055464 0.0158482 -0.350 0.727
M5 -0.0061031 0.0158215 -0.386 0.701
M6 -0.0091598 0.0157983 -0.580 0.564
M7 -0.0097165 0.0157787 -0.616 0.540
M8 -0.0052732 0.0157626 -0.335 0.739
M9 -0.0020799 0.0157501 -0.132 0.895
M10 0.0023634 0.0157411 0.150 0.881
M11 0.0005567 0.0157357 0.035 0.972
t 0.0030567 0.0002374 12.875 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.03147 on 82 degrees of freedom
Multiple R-squared: 0.939, Adjusted R-squared: 0.9294
F-statistic: 97.18 on 13 and 82 DF, p-value: < 2.2e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 4.925269e-43 9.850538e-43 1.0000000
[2,] 1.314742e-56 2.629485e-56 1.0000000
[3,] 2.422119e-05 4.844239e-05 0.9999758
[4,] 2.422993e-02 4.845986e-02 0.9757701
[5,] 5.403111e-02 1.080622e-01 0.9459689
[6,] 7.188277e-02 1.437655e-01 0.9281172
[7,] 8.154520e-02 1.630904e-01 0.9184548
[8,] 8.583876e-02 1.716775e-01 0.9141612
[9,] 8.155843e-02 1.631169e-01 0.9184416
[10,] 6.873316e-02 1.374663e-01 0.9312668
[11,] 5.723581e-02 1.144716e-01 0.9427642
[12,] 4.743619e-02 9.487238e-02 0.9525638
[13,] 3.926722e-02 7.853444e-02 0.9607328
[14,] 3.424423e-02 6.848846e-02 0.9657558
[15,] 3.021534e-02 6.043067e-02 0.9697847
[16,] 3.496209e-02 6.992418e-02 0.9650379
[17,] 3.642053e-02 7.284106e-02 0.9635795
[18,] 3.327432e-02 6.654864e-02 0.9667257
[19,] 3.021834e-02 6.043669e-02 0.9697817
[20,] 2.668417e-02 5.336834e-02 0.9733158
[21,] 3.220574e-02 6.441149e-02 0.9677943
[22,] 2.969643e-02 5.939286e-02 0.9703036
[23,] 2.595176e-02 5.190353e-02 0.9740482
[24,] 2.158528e-02 4.317055e-02 0.9784147
[25,] 1.712052e-02 3.424103e-02 0.9828795
[26,] 1.387001e-02 2.774001e-02 0.9861300
[27,] 1.206095e-02 2.412189e-02 0.9879391
[28,] 1.547769e-02 3.095538e-02 0.9845223
[29,] 1.640420e-02 3.280840e-02 0.9835958
[30,] 1.520650e-02 3.041300e-02 0.9847935
[31,] 1.337274e-02 2.674548e-02 0.9866273
[32,] 1.154576e-02 2.309152e-02 0.9884542
[33,] 9.587537e-03 1.917507e-02 0.9904125
[34,] 6.404849e-03 1.280970e-02 0.9935952
[35,] 4.904641e-03 9.809283e-03 0.9950954
[36,] 3.329192e-03 6.658384e-03 0.9966708
[37,] 2.060839e-03 4.121678e-03 0.9979392
[38,] 1.253952e-03 2.507905e-03 0.9987460
[39,] 9.688969e-04 1.937794e-03 0.9990311
[40,] 9.392664e-04 1.878533e-03 0.9990607
[41,] 6.264166e-04 1.252833e-03 0.9993736
[42,] 4.578530e-04 9.157061e-04 0.9995421
[43,] 3.063202e-04 6.126405e-04 0.9996937
[44,] 1.854239e-04 3.708477e-04 0.9998146
[45,] 2.664785e-04 5.329569e-04 0.9997335
[46,] 6.062150e-04 1.212430e-03 0.9993938
[47,] 1.021701e-03 2.043402e-03 0.9989783
[48,] 3.343228e-03 6.686456e-03 0.9966568
[49,] 6.933520e-03 1.386704e-02 0.9930665
[50,] 1.485121e-02 2.970242e-02 0.9851488
[51,] 3.289871e-02 6.579742e-02 0.9671013
[52,] 6.343835e-02 1.268767e-01 0.9365617
[53,] 1.113536e-01 2.227072e-01 0.8886464
[54,] 9.850264e-02 1.970053e-01 0.9014974
[55,] 8.108858e-02 1.621772e-01 0.9189114
[56,] 5.815499e-02 1.163100e-01 0.9418450
[57,] 4.855814e-02 9.711628e-02 0.9514419
[58,] 4.566157e-02 9.132314e-02 0.9543384
[59,] 3.363459e-02 6.726918e-02 0.9663654
[60,] 3.641819e-02 7.283638e-02 0.9635818
[61,] 4.956327e-02 9.912654e-02 0.9504367
[62,] 7.815352e-02 1.563070e-01 0.9218465
[63,] 1.713672e-01 3.427343e-01 0.8286328
> postscript(file="/var/www/html/freestat/rcomp/tmp/1ml1b1227816104.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/2j9m71227816104.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/3pqi81227816104.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/4ssjz1227816104.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/5jjmt1227816104.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 = 96
Frequency = 1
1 2 3 4 5
0.0463711269 0.0394511201 0.0369511201 0.0344511201 0.0319511201
6 7 8 9 10
0.0319511201 0.0294511201 0.0219511201 0.0157011201 0.0082011201
11 12 13 14 15
0.0069511201 0.0044511201 0.0096907991 0.0027707922 0.0002707922
16 17 18 19 20
-0.0022292078 -0.0047292078 -0.0047292078 0.0027707922 0.0352707922
21 22 23 24 25
0.0290207922 0.0215207922 0.0202707922 0.0177707922 0.0230104712
26 27 28 29 30
0.0160904643 0.0135904643 0.0110904643 0.0085904643 0.0085904643
31 32 33 34 35
0.0060904643 -0.0014095357 -0.0076595357 -0.0151595357 -0.0164095357
36 37 38 39 40
-0.0189095357 -0.0136698567 -0.0205898635 -0.0230898635 -0.0255898635
41 42 43 44 45
-0.0280898635 -0.0280898635 -0.0305898635 -0.0380898635 -0.0443398635
46 47 48 49 50
-0.0518398635 -0.0530898635 -0.0555898635 -0.0503501845 -0.0319101365
51 52 53 54 55
-0.0244101365 -0.0269101365 -0.0294101365 -0.0294101365 -0.0219101365
56 57 58 59 60
-0.0194101365 -0.0256601365 -0.0231601365 -0.0244101365 -0.0269101365
61 62 63 64 65
-0.0116704575 -0.0085904643 -0.0110904643 -0.0035904643 -0.0060904643
66 67 68 69 70
-0.0060904643 -0.0085904643 -0.0160904643 -0.0123404643 -0.0198404643
71 72 73 74 75
-0.0210904643 -0.0235904643 -0.0183507853 -0.0152707922 -0.0177707922
76 77 78 79 80
-0.0102707922 -0.0027707922 -0.0027707922 -0.0052707922 -0.0127707922
81 82 83 84 85
-0.0090207922 -0.0065207922 -0.0077707922 -0.0002707922 0.0149688868
86 87 88 89 90
0.0180488799 0.0255488799 0.0230488799 0.0305488799 0.0305488799
91 92 93 94 95
0.0280488799 0.0305488799 0.0542988799 0.0867988799 0.0955488799
96
0.1030488799
> postscript(file="/var/www/html/freestat/rcomp/tmp/609rr1227816104.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 = 96
Frequency = 1
lag(myerror, k = 1) myerror
0 0.0463711269 NA
1 0.0394511201 0.0463711269
2 0.0369511201 0.0394511201
3 0.0344511201 0.0369511201
4 0.0319511201 0.0344511201
5 0.0319511201 0.0319511201
6 0.0294511201 0.0319511201
7 0.0219511201 0.0294511201
8 0.0157011201 0.0219511201
9 0.0082011201 0.0157011201
10 0.0069511201 0.0082011201
11 0.0044511201 0.0069511201
12 0.0096907991 0.0044511201
13 0.0027707922 0.0096907991
14 0.0002707922 0.0027707922
15 -0.0022292078 0.0002707922
16 -0.0047292078 -0.0022292078
17 -0.0047292078 -0.0047292078
18 0.0027707922 -0.0047292078
19 0.0352707922 0.0027707922
20 0.0290207922 0.0352707922
21 0.0215207922 0.0290207922
22 0.0202707922 0.0215207922
23 0.0177707922 0.0202707922
24 0.0230104712 0.0177707922
25 0.0160904643 0.0230104712
26 0.0135904643 0.0160904643
27 0.0110904643 0.0135904643
28 0.0085904643 0.0110904643
29 0.0085904643 0.0085904643
30 0.0060904643 0.0085904643
31 -0.0014095357 0.0060904643
32 -0.0076595357 -0.0014095357
33 -0.0151595357 -0.0076595357
34 -0.0164095357 -0.0151595357
35 -0.0189095357 -0.0164095357
36 -0.0136698567 -0.0189095357
37 -0.0205898635 -0.0136698567
38 -0.0230898635 -0.0205898635
39 -0.0255898635 -0.0230898635
40 -0.0280898635 -0.0255898635
41 -0.0280898635 -0.0280898635
42 -0.0305898635 -0.0280898635
43 -0.0380898635 -0.0305898635
44 -0.0443398635 -0.0380898635
45 -0.0518398635 -0.0443398635
46 -0.0530898635 -0.0518398635
47 -0.0555898635 -0.0530898635
48 -0.0503501845 -0.0555898635
49 -0.0319101365 -0.0503501845
50 -0.0244101365 -0.0319101365
51 -0.0269101365 -0.0244101365
52 -0.0294101365 -0.0269101365
53 -0.0294101365 -0.0294101365
54 -0.0219101365 -0.0294101365
55 -0.0194101365 -0.0219101365
56 -0.0256601365 -0.0194101365
57 -0.0231601365 -0.0256601365
58 -0.0244101365 -0.0231601365
59 -0.0269101365 -0.0244101365
60 -0.0116704575 -0.0269101365
61 -0.0085904643 -0.0116704575
62 -0.0110904643 -0.0085904643
63 -0.0035904643 -0.0110904643
64 -0.0060904643 -0.0035904643
65 -0.0060904643 -0.0060904643
66 -0.0085904643 -0.0060904643
67 -0.0160904643 -0.0085904643
68 -0.0123404643 -0.0160904643
69 -0.0198404643 -0.0123404643
70 -0.0210904643 -0.0198404643
71 -0.0235904643 -0.0210904643
72 -0.0183507853 -0.0235904643
73 -0.0152707922 -0.0183507853
74 -0.0177707922 -0.0152707922
75 -0.0102707922 -0.0177707922
76 -0.0027707922 -0.0102707922
77 -0.0027707922 -0.0027707922
78 -0.0052707922 -0.0027707922
79 -0.0127707922 -0.0052707922
80 -0.0090207922 -0.0127707922
81 -0.0065207922 -0.0090207922
82 -0.0077707922 -0.0065207922
83 -0.0002707922 -0.0077707922
84 0.0149688868 -0.0002707922
85 0.0180488799 0.0149688868
86 0.0255488799 0.0180488799
87 0.0230488799 0.0255488799
88 0.0305488799 0.0230488799
89 0.0305488799 0.0305488799
90 0.0280488799 0.0305488799
91 0.0305488799 0.0280488799
92 0.0542988799 0.0305488799
93 0.0867988799 0.0542988799
94 0.0955488799 0.0867988799
95 0.1030488799 0.0955488799
96 NA 0.1030488799
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.0394511201 0.0463711269
[2,] 0.0369511201 0.0394511201
[3,] 0.0344511201 0.0369511201
[4,] 0.0319511201 0.0344511201
[5,] 0.0319511201 0.0319511201
[6,] 0.0294511201 0.0319511201
[7,] 0.0219511201 0.0294511201
[8,] 0.0157011201 0.0219511201
[9,] 0.0082011201 0.0157011201
[10,] 0.0069511201 0.0082011201
[11,] 0.0044511201 0.0069511201
[12,] 0.0096907991 0.0044511201
[13,] 0.0027707922 0.0096907991
[14,] 0.0002707922 0.0027707922
[15,] -0.0022292078 0.0002707922
[16,] -0.0047292078 -0.0022292078
[17,] -0.0047292078 -0.0047292078
[18,] 0.0027707922 -0.0047292078
[19,] 0.0352707922 0.0027707922
[20,] 0.0290207922 0.0352707922
[21,] 0.0215207922 0.0290207922
[22,] 0.0202707922 0.0215207922
[23,] 0.0177707922 0.0202707922
[24,] 0.0230104712 0.0177707922
[25,] 0.0160904643 0.0230104712
[26,] 0.0135904643 0.0160904643
[27,] 0.0110904643 0.0135904643
[28,] 0.0085904643 0.0110904643
[29,] 0.0085904643 0.0085904643
[30,] 0.0060904643 0.0085904643
[31,] -0.0014095357 0.0060904643
[32,] -0.0076595357 -0.0014095357
[33,] -0.0151595357 -0.0076595357
[34,] -0.0164095357 -0.0151595357
[35,] -0.0189095357 -0.0164095357
[36,] -0.0136698567 -0.0189095357
[37,] -0.0205898635 -0.0136698567
[38,] -0.0230898635 -0.0205898635
[39,] -0.0255898635 -0.0230898635
[40,] -0.0280898635 -0.0255898635
[41,] -0.0280898635 -0.0280898635
[42,] -0.0305898635 -0.0280898635
[43,] -0.0380898635 -0.0305898635
[44,] -0.0443398635 -0.0380898635
[45,] -0.0518398635 -0.0443398635
[46,] -0.0530898635 -0.0518398635
[47,] -0.0555898635 -0.0530898635
[48,] -0.0503501845 -0.0555898635
[49,] -0.0319101365 -0.0503501845
[50,] -0.0244101365 -0.0319101365
[51,] -0.0269101365 -0.0244101365
[52,] -0.0294101365 -0.0269101365
[53,] -0.0294101365 -0.0294101365
[54,] -0.0219101365 -0.0294101365
[55,] -0.0194101365 -0.0219101365
[56,] -0.0256601365 -0.0194101365
[57,] -0.0231601365 -0.0256601365
[58,] -0.0244101365 -0.0231601365
[59,] -0.0269101365 -0.0244101365
[60,] -0.0116704575 -0.0269101365
[61,] -0.0085904643 -0.0116704575
[62,] -0.0110904643 -0.0085904643
[63,] -0.0035904643 -0.0110904643
[64,] -0.0060904643 -0.0035904643
[65,] -0.0060904643 -0.0060904643
[66,] -0.0085904643 -0.0060904643
[67,] -0.0160904643 -0.0085904643
[68,] -0.0123404643 -0.0160904643
[69,] -0.0198404643 -0.0123404643
[70,] -0.0210904643 -0.0198404643
[71,] -0.0235904643 -0.0210904643
[72,] -0.0183507853 -0.0235904643
[73,] -0.0152707922 -0.0183507853
[74,] -0.0177707922 -0.0152707922
[75,] -0.0102707922 -0.0177707922
[76,] -0.0027707922 -0.0102707922
[77,] -0.0027707922 -0.0027707922
[78,] -0.0052707922 -0.0027707922
[79,] -0.0127707922 -0.0052707922
[80,] -0.0090207922 -0.0127707922
[81,] -0.0065207922 -0.0090207922
[82,] -0.0077707922 -0.0065207922
[83,] -0.0002707922 -0.0077707922
[84,] 0.0149688868 -0.0002707922
[85,] 0.0180488799 0.0149688868
[86,] 0.0255488799 0.0180488799
[87,] 0.0230488799 0.0255488799
[88,] 0.0305488799 0.0230488799
[89,] 0.0305488799 0.0305488799
[90,] 0.0280488799 0.0305488799
[91,] 0.0305488799 0.0280488799
[92,] 0.0542988799 0.0305488799
[93,] 0.0867988799 0.0542988799
[94,] 0.0955488799 0.0867988799
[95,] 0.1030488799 0.0955488799
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.0394511201 0.0463711269
2 0.0369511201 0.0394511201
3 0.0344511201 0.0369511201
4 0.0319511201 0.0344511201
5 0.0319511201 0.0319511201
6 0.0294511201 0.0319511201
7 0.0219511201 0.0294511201
8 0.0157011201 0.0219511201
9 0.0082011201 0.0157011201
10 0.0069511201 0.0082011201
11 0.0044511201 0.0069511201
12 0.0096907991 0.0044511201
13 0.0027707922 0.0096907991
14 0.0002707922 0.0027707922
15 -0.0022292078 0.0002707922
16 -0.0047292078 -0.0022292078
17 -0.0047292078 -0.0047292078
18 0.0027707922 -0.0047292078
19 0.0352707922 0.0027707922
20 0.0290207922 0.0352707922
21 0.0215207922 0.0290207922
22 0.0202707922 0.0215207922
23 0.0177707922 0.0202707922
24 0.0230104712 0.0177707922
25 0.0160904643 0.0230104712
26 0.0135904643 0.0160904643
27 0.0110904643 0.0135904643
28 0.0085904643 0.0110904643
29 0.0085904643 0.0085904643
30 0.0060904643 0.0085904643
31 -0.0014095357 0.0060904643
32 -0.0076595357 -0.0014095357
33 -0.0151595357 -0.0076595357
34 -0.0164095357 -0.0151595357
35 -0.0189095357 -0.0164095357
36 -0.0136698567 -0.0189095357
37 -0.0205898635 -0.0136698567
38 -0.0230898635 -0.0205898635
39 -0.0255898635 -0.0230898635
40 -0.0280898635 -0.0255898635
41 -0.0280898635 -0.0280898635
42 -0.0305898635 -0.0280898635
43 -0.0380898635 -0.0305898635
44 -0.0443398635 -0.0380898635
45 -0.0518398635 -0.0443398635
46 -0.0530898635 -0.0518398635
47 -0.0555898635 -0.0530898635
48 -0.0503501845 -0.0555898635
49 -0.0319101365 -0.0503501845
50 -0.0244101365 -0.0319101365
51 -0.0269101365 -0.0244101365
52 -0.0294101365 -0.0269101365
53 -0.0294101365 -0.0294101365
54 -0.0219101365 -0.0294101365
55 -0.0194101365 -0.0219101365
56 -0.0256601365 -0.0194101365
57 -0.0231601365 -0.0256601365
58 -0.0244101365 -0.0231601365
59 -0.0269101365 -0.0244101365
60 -0.0116704575 -0.0269101365
61 -0.0085904643 -0.0116704575
62 -0.0110904643 -0.0085904643
63 -0.0035904643 -0.0110904643
64 -0.0060904643 -0.0035904643
65 -0.0060904643 -0.0060904643
66 -0.0085904643 -0.0060904643
67 -0.0160904643 -0.0085904643
68 -0.0123404643 -0.0160904643
69 -0.0198404643 -0.0123404643
70 -0.0210904643 -0.0198404643
71 -0.0235904643 -0.0210904643
72 -0.0183507853 -0.0235904643
73 -0.0152707922 -0.0183507853
74 -0.0177707922 -0.0152707922
75 -0.0102707922 -0.0177707922
76 -0.0027707922 -0.0102707922
77 -0.0027707922 -0.0027707922
78 -0.0052707922 -0.0027707922
79 -0.0127707922 -0.0052707922
80 -0.0090207922 -0.0127707922
81 -0.0065207922 -0.0090207922
82 -0.0077707922 -0.0065207922
83 -0.0002707922 -0.0077707922
84 0.0149688868 -0.0002707922
85 0.0180488799 0.0149688868
86 0.0255488799 0.0180488799
87 0.0230488799 0.0255488799
88 0.0305488799 0.0230488799
89 0.0305488799 0.0305488799
90 0.0280488799 0.0305488799
91 0.0305488799 0.0280488799
92 0.0542988799 0.0305488799
93 0.0867988799 0.0542988799
94 0.0955488799 0.0867988799
95 0.1030488799 0.0955488799
> 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/7wr6f1227816104.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/8twvc1227816104.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/9i7za1227816104.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/10mkhu1227816104.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/11zh161227816104.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/12j0po1227816104.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/13m7yf1227816104.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/141lvd1227816104.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/15s2lg1227816104.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/16oxld1227816105.tab")
+ }
>
> system("convert tmp/1ml1b1227816104.ps tmp/1ml1b1227816104.png")
> system("convert tmp/2j9m71227816104.ps tmp/2j9m71227816104.png")
> system("convert tmp/3pqi81227816104.ps tmp/3pqi81227816104.png")
> system("convert tmp/4ssjz1227816104.ps tmp/4ssjz1227816104.png")
> system("convert tmp/5jjmt1227816104.ps tmp/5jjmt1227816104.png")
> system("convert tmp/609rr1227816104.ps tmp/609rr1227816104.png")
> system("convert tmp/7wr6f1227816104.ps tmp/7wr6f1227816104.png")
> system("convert tmp/8twvc1227816104.ps tmp/8twvc1227816104.png")
> system("convert tmp/9i7za1227816104.ps tmp/9i7za1227816104.png")
> system("convert tmp/10mkhu1227816104.ps tmp/10mkhu1227816104.png")
>
>
> proc.time()
user system elapsed
4.160 2.493 4.493