R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(98.1
+ ,102.8
+ ,104.7
+ ,95.9
+ ,94.6
+ ,113.9
+ ,98.1
+ ,102.8
+ ,104.7
+ ,95.9
+ ,80.9
+ ,113.9
+ ,98.1
+ ,102.8
+ ,104.7
+ ,95.7
+ ,80.9
+ ,113.9
+ ,98.1
+ ,102.8
+ ,113.2
+ ,95.7
+ ,80.9
+ ,113.9
+ ,98.1
+ ,105.9
+ ,113.2
+ ,95.7
+ ,80.9
+ ,113.9
+ ,108.8
+ ,105.9
+ ,113.2
+ ,95.7
+ ,80.9
+ ,102.3
+ ,108.8
+ ,105.9
+ ,113.2
+ ,95.7
+ ,99
+ ,102.3
+ ,108.8
+ ,105.9
+ ,113.2
+ ,100.7
+ ,99
+ ,102.3
+ ,108.8
+ ,105.9
+ ,115.5
+ ,100.7
+ ,99
+ ,102.3
+ ,108.8
+ ,100.7
+ ,115.5
+ ,100.7
+ ,99
+ ,102.3
+ ,109.9
+ ,100.7
+ ,115.5
+ ,100.7
+ ,99
+ ,114.6
+ ,109.9
+ ,100.7
+ ,115.5
+ ,100.7
+ ,85.4
+ ,114.6
+ ,109.9
+ ,100.7
+ ,115.5
+ ,100.5
+ ,85.4
+ ,114.6
+ ,109.9
+ ,100.7
+ ,114.8
+ ,100.5
+ ,85.4
+ ,114.6
+ ,109.9
+ ,116.5
+ ,114.8
+ ,100.5
+ ,85.4
+ ,114.6
+ ,112.9
+ ,116.5
+ ,114.8
+ ,100.5
+ ,85.4
+ ,102
+ ,112.9
+ ,116.5
+ ,114.8
+ ,100.5
+ ,106
+ ,102
+ ,112.9
+ ,116.5
+ ,114.8
+ ,105.3
+ ,106
+ ,102
+ ,112.9
+ ,116.5
+ ,118.8
+ ,105.3
+ ,106
+ ,102
+ ,112.9
+ ,106.1
+ ,118.8
+ ,105.3
+ ,106
+ ,102
+ ,109.3
+ ,106.1
+ ,118.8
+ ,105.3
+ ,106
+ ,117.2
+ ,109.3
+ ,106.1
+ ,118.8
+ ,105.3
+ ,92.5
+ ,117.2
+ ,109.3
+ ,106.1
+ ,118.8
+ ,104.2
+ ,92.5
+ ,117.2
+ ,109.3
+ ,106.1
+ ,112.5
+ ,104.2
+ ,92.5
+ ,117.2
+ ,109.3
+ ,122.4
+ ,112.5
+ ,104.2
+ ,92.5
+ ,117.2
+ ,113.3
+ ,122.4
+ ,112.5
+ ,104.2
+ ,92.5
+ ,100
+ ,113.3
+ ,122.4
+ ,112.5
+ ,104.2
+ ,110.7
+ ,100
+ ,113.3
+ ,122.4
+ ,112.5
+ ,112.8
+ ,110.7
+ ,100
+ ,113.3
+ ,122.4
+ ,109.8
+ ,112.8
+ ,110.7
+ ,100
+ ,113.3
+ ,117.3
+ ,109.8
+ ,112.8
+ ,110.7
+ ,100
+ ,109.1
+ ,117.3
+ ,109.8
+ ,112.8
+ ,110.7
+ ,115.9
+ ,109.1
+ ,117.3
+ ,109.8
+ ,112.8
+ ,96
+ ,115.9
+ ,109.1
+ ,117.3
+ ,109.8
+ ,99.8
+ ,96
+ ,115.9
+ ,109.1
+ ,117.3
+ ,116.8
+ ,99.8
+ ,96
+ ,115.9
+ ,109.1
+ ,115.7
+ ,116.8
+ ,99.8
+ ,96
+ ,115.9
+ ,99.4
+ ,115.7
+ ,116.8
+ ,99.8
+ ,96
+ ,94.3
+ ,99.4
+ ,115.7
+ ,116.8
+ ,99.8
+ ,91
+ ,94.3
+ ,99.4
+ ,115.7
+ ,116.8
+ ,93.2
+ ,91
+ ,94.3
+ ,99.4
+ ,115.7
+ ,103.1
+ ,93.2
+ ,91
+ ,94.3
+ ,99.4
+ ,94.1
+ ,103.1
+ ,93.2
+ ,91
+ ,94.3
+ ,91.8
+ ,94.1
+ ,103.1
+ ,93.2
+ ,91
+ ,102.7
+ ,91.8
+ ,94.1
+ ,103.1
+ ,93.2
+ ,82.6
+ ,102.7
+ ,91.8
+ ,94.1
+ ,103.1
+ ,89.1
+ ,82.6
+ ,102.7
+ ,91.8
+ ,94.1
+ ,104.5
+ ,89.1
+ ,82.6
+ ,102.7
+ ,91.8
+ ,105.1
+ ,104.5
+ ,89.1
+ ,82.6
+ ,102.7
+ ,95.1
+ ,105.1
+ ,104.5
+ ,89.1
+ ,82.6
+ ,88.7
+ ,95.1
+ ,105.1
+ ,104.5
+ ,89.1
+ ,86.3
+ ,88.7
+ ,95.1
+ ,105.1
+ ,104.5
+ ,91.8
+ ,86.3
+ ,88.7
+ ,95.1
+ ,105.1
+ ,111.5
+ ,91.8
+ ,86.3
+ ,88.7
+ ,95.1
+ ,99.7
+ ,111.5
+ ,91.8
+ ,86.3
+ ,88.7
+ ,97.5
+ ,99.7
+ ,111.5
+ ,91.8
+ ,86.3
+ ,111.7
+ ,97.5
+ ,99.7
+ ,111.5
+ ,91.8
+ ,86.2
+ ,111.7
+ ,97.5
+ ,99.7
+ ,111.5
+ ,95.4
+ ,86.2
+ ,111.7
+ ,97.5
+ ,99.7)
+ ,dim=c(5
+ ,64)
+ ,dimnames=list(c('y'
+ ,'y1'
+ ,'y2'
+ ,'y3'
+ ,'y4')
+ ,1:64))
> y <- array(NA,dim=c(5,64),dimnames=list(c('y','y1','y2','y3','y4'),1:64))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = '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 y1 y2 y3 y4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 98.1 102.8 104.7 95.9 94.6 1 0 0 0 0 0 0 0 0 0 0 1
2 113.9 98.1 102.8 104.7 95.9 0 1 0 0 0 0 0 0 0 0 0 2
3 80.9 113.9 98.1 102.8 104.7 0 0 1 0 0 0 0 0 0 0 0 3
4 95.7 80.9 113.9 98.1 102.8 0 0 0 1 0 0 0 0 0 0 0 4
5 113.2 95.7 80.9 113.9 98.1 0 0 0 0 1 0 0 0 0 0 0 5
6 105.9 113.2 95.7 80.9 113.9 0 0 0 0 0 1 0 0 0 0 0 6
7 108.8 105.9 113.2 95.7 80.9 0 0 0 0 0 0 1 0 0 0 0 7
8 102.3 108.8 105.9 113.2 95.7 0 0 0 0 0 0 0 1 0 0 0 8
9 99.0 102.3 108.8 105.9 113.2 0 0 0 0 0 0 0 0 1 0 0 9
10 100.7 99.0 102.3 108.8 105.9 0 0 0 0 0 0 0 0 0 1 0 10
11 115.5 100.7 99.0 102.3 108.8 0 0 0 0 0 0 0 0 0 0 1 11
12 100.7 115.5 100.7 99.0 102.3 0 0 0 0 0 0 0 0 0 0 0 12
13 109.9 100.7 115.5 100.7 99.0 1 0 0 0 0 0 0 0 0 0 0 13
14 114.6 109.9 100.7 115.5 100.7 0 1 0 0 0 0 0 0 0 0 0 14
15 85.4 114.6 109.9 100.7 115.5 0 0 1 0 0 0 0 0 0 0 0 15
16 100.5 85.4 114.6 109.9 100.7 0 0 0 1 0 0 0 0 0 0 0 16
17 114.8 100.5 85.4 114.6 109.9 0 0 0 0 1 0 0 0 0 0 0 17
18 116.5 114.8 100.5 85.4 114.6 0 0 0 0 0 1 0 0 0 0 0 18
19 112.9 116.5 114.8 100.5 85.4 0 0 0 0 0 0 1 0 0 0 0 19
20 102.0 112.9 116.5 114.8 100.5 0 0 0 0 0 0 0 1 0 0 0 20
21 106.0 102.0 112.9 116.5 114.8 0 0 0 0 0 0 0 0 1 0 0 21
22 105.3 106.0 102.0 112.9 116.5 0 0 0 0 0 0 0 0 0 1 0 22
23 118.8 105.3 106.0 102.0 112.9 0 0 0 0 0 0 0 0 0 0 1 23
24 106.1 118.8 105.3 106.0 102.0 0 0 0 0 0 0 0 0 0 0 0 24
25 109.3 106.1 118.8 105.3 106.0 1 0 0 0 0 0 0 0 0 0 0 25
26 117.2 109.3 106.1 118.8 105.3 0 1 0 0 0 0 0 0 0 0 0 26
27 92.5 117.2 109.3 106.1 118.8 0 0 1 0 0 0 0 0 0 0 0 27
28 104.2 92.5 117.2 109.3 106.1 0 0 0 1 0 0 0 0 0 0 0 28
29 112.5 104.2 92.5 117.2 109.3 0 0 0 0 1 0 0 0 0 0 0 29
30 122.4 112.5 104.2 92.5 117.2 0 0 0 0 0 1 0 0 0 0 0 30
31 113.3 122.4 112.5 104.2 92.5 0 0 0 0 0 0 1 0 0 0 0 31
32 100.0 113.3 122.4 112.5 104.2 0 0 0 0 0 0 0 1 0 0 0 32
33 110.7 100.0 113.3 122.4 112.5 0 0 0 0 0 0 0 0 1 0 0 33
34 112.8 110.7 100.0 113.3 122.4 0 0 0 0 0 0 0 0 0 1 0 34
35 109.8 112.8 110.7 100.0 113.3 0 0 0 0 0 0 0 0 0 0 1 35
36 117.3 109.8 112.8 110.7 100.0 0 0 0 0 0 0 0 0 0 0 0 36
37 109.1 117.3 109.8 112.8 110.7 1 0 0 0 0 0 0 0 0 0 0 37
38 115.9 109.1 117.3 109.8 112.8 0 1 0 0 0 0 0 0 0 0 0 38
39 96.0 115.9 109.1 117.3 109.8 0 0 1 0 0 0 0 0 0 0 0 39
40 99.8 96.0 115.9 109.1 117.3 0 0 0 1 0 0 0 0 0 0 0 40
41 116.8 99.8 96.0 115.9 109.1 0 0 0 0 1 0 0 0 0 0 0 41
42 115.7 116.8 99.8 96.0 115.9 0 0 0 0 0 1 0 0 0 0 0 42
43 99.4 115.7 116.8 99.8 96.0 0 0 0 0 0 0 1 0 0 0 0 43
44 94.3 99.4 115.7 116.8 99.8 0 0 0 0 0 0 0 1 0 0 0 44
45 91.0 94.3 99.4 115.7 116.8 0 0 0 0 0 0 0 0 1 0 0 45
46 93.2 91.0 94.3 99.4 115.7 0 0 0 0 0 0 0 0 0 1 0 46
47 103.1 93.2 91.0 94.3 99.4 0 0 0 0 0 0 0 0 0 0 1 47
48 94.1 103.1 93.2 91.0 94.3 0 0 0 0 0 0 0 0 0 0 0 48
49 91.8 94.1 103.1 93.2 91.0 1 0 0 0 0 0 0 0 0 0 0 49
50 102.7 91.8 94.1 103.1 93.2 0 1 0 0 0 0 0 0 0 0 0 50
51 82.6 102.7 91.8 94.1 103.1 0 0 1 0 0 0 0 0 0 0 0 51
52 89.1 82.6 102.7 91.8 94.1 0 0 0 1 0 0 0 0 0 0 0 52
53 104.5 89.1 82.6 102.7 91.8 0 0 0 0 1 0 0 0 0 0 0 53
54 105.1 104.5 89.1 82.6 102.7 0 0 0 0 0 1 0 0 0 0 0 54
55 95.1 105.1 104.5 89.1 82.6 0 0 0 0 0 0 1 0 0 0 0 55
56 88.7 95.1 105.1 104.5 89.1 0 0 0 0 0 0 0 1 0 0 0 56
57 86.3 88.7 95.1 105.1 104.5 0 0 0 0 0 0 0 0 1 0 0 57
58 91.8 86.3 88.7 95.1 105.1 0 0 0 0 0 0 0 0 0 1 0 58
59 111.5 91.8 86.3 88.7 95.1 0 0 0 0 0 0 0 0 0 0 1 59
60 99.7 111.5 91.8 86.3 88.7 0 0 0 0 0 0 0 0 0 0 0 60
61 97.5 99.7 111.5 91.8 86.3 1 0 0 0 0 0 0 0 0 0 0 61
62 111.7 97.5 99.7 111.5 91.8 0 1 0 0 0 0 0 0 0 0 0 62
63 86.2 111.7 97.5 99.7 111.5 0 0 1 0 0 0 0 0 0 0 0 63
64 95.4 86.2 111.7 97.5 99.7 0 0 0 1 0 0 0 0 0 0 0 64
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) y1 y2 y3 y4 M1
18.18644 0.20675 0.33378 0.52082 -0.22036 -3.29046
M2 M3 M4 M5 M6 M7
4.25601 -16.86385 -5.88390 9.78051 19.36206 -3.81998
M8 M9 M10 M11 t
-16.25169 -8.12360 0.57942 13.57088 -0.03381
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-8.16021 -2.25177 0.09915 2.53637 8.71187
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 18.18644 11.16405 1.629 0.109995
y1 0.20675 0.14245 1.451 0.153298
y2 0.33378 0.12221 2.731 0.008859 **
y3 0.52082 0.12073 4.314 8.18e-05 ***
y4 -0.22036 0.14010 -1.573 0.122458
M1 -3.29046 3.05909 -1.076 0.287582
M2 4.25601 3.37981 1.259 0.214160
M3 -16.86385 2.99816 -5.625 9.91e-07 ***
M4 -5.88390 5.00468 -1.176 0.245644
M5 9.78051 4.76462 2.053 0.045689 *
M6 19.36206 3.80041 5.095 6.10e-06 ***
M7 -3.81998 3.50489 -1.090 0.281313
M8 -16.25169 3.46892 -4.685 2.42e-05 ***
M9 -8.12360 4.78888 -1.696 0.096435 .
M10 0.57942 4.62382 0.125 0.900811
M11 13.57088 3.55086 3.822 0.000388 ***
t -0.03381 0.03105 -1.089 0.281733
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.099 on 47 degrees of freedom
Multiple R-squared: 0.8785, Adjusted R-squared: 0.8372
F-statistic: 21.24 on 16 and 47 DF, p-value: 2.985e-16
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 0.285184637 0.57036927 0.7148154
[2,] 0.168978882 0.33795776 0.8310211
[3,] 0.091458957 0.18291791 0.9085410
[4,] 0.051465475 0.10293095 0.9485345
[5,] 0.030353957 0.06070791 0.9696460
[6,] 0.020733439 0.04146688 0.9792666
[7,] 0.035995190 0.07199038 0.9640048
[8,] 0.017587839 0.03517568 0.9824122
[9,] 0.007904095 0.01580819 0.9920959
[10,] 0.032316457 0.06463291 0.9676835
[11,] 0.036380123 0.07276025 0.9636199
[12,] 0.025782896 0.05156579 0.9742171
[13,] 0.043839219 0.08767844 0.9561608
[14,] 0.088284109 0.17656822 0.9117159
[15,] 0.160033594 0.32006719 0.8399664
[16,] 0.269534321 0.53906864 0.7304657
[17,] 0.572302033 0.85539593 0.4276980
[18,] 0.517358034 0.96528393 0.4826420
[19,] 0.464070351 0.92814070 0.5359296
[20,] 0.363251568 0.72650314 0.6367484
[21,] 0.315468915 0.63093783 0.6845311
[22,] 0.336845561 0.67369112 0.6631544
[23,] 0.344900622 0.68980124 0.6550994
[24,] 0.454425030 0.90885006 0.5455750
[25,] 0.391720319 0.78344064 0.6082797
> postscript(file="/var/www/html/freestat/rcomp/tmp/1j90y1291055503.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/2j90y1291055503.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/3j90y1291055503.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/4t0hj1291055503.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/5t0hj1291055503.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 = 64
Frequency = 1
1 2 3 4 5 6 7
-2.0635621 3.5329087 -7.0825968 0.3495696 0.9089827 -3.8279076 2.9759768
8 9 10 11 12 13 14
4.9254253 1.5654259 -5.6709834 0.9457597 -3.5905325 5.4411523 -1.6672936
15 16 17 18 19 20 21
-2.7865613 -2.2172448 2.6559675 3.0553994 3.2477355 0.8698169 2.4965039
22 23 24 25 26 27 28
-1.8119560 2.4236875 -3.7143702 2.1756675 -1.0449988 2.2965927 1.0551433
29 30 31 32 33 34 35
-3.8594973 5.4767383 3.2388077 -0.7632393 4.3024987 6.8813701 -8.1602143
36 37 38 39 40 41 42
4.3601626 0.1992655 0.7038718 -1.2786653 -0.6566598 1.2206849 -2.3473297
43 44 45 46 47 48 49
-7.2426091 -4.1565198 -4.7367314 -0.5743614 -3.9210982 -1.5026892 -3.7950658
50 51 52 53 54 55 56
-1.5995584 4.8372026 0.1232347 -0.9261378 -2.3569004 -2.2199110 -0.8754832
57 58 59 60 61 62 63
-3.6276971 1.1759307 8.7118653 4.4474293 -1.9574574 0.0750702 4.0140281
64
1.3459570
> postscript(file="/var/www/html/freestat/rcomp/tmp/6mszm1291055503.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 = 64
Frequency = 1
lag(myerror, k = 1) myerror
0 -2.0635621 NA
1 3.5329087 -2.0635621
2 -7.0825968 3.5329087
3 0.3495696 -7.0825968
4 0.9089827 0.3495696
5 -3.8279076 0.9089827
6 2.9759768 -3.8279076
7 4.9254253 2.9759768
8 1.5654259 4.9254253
9 -5.6709834 1.5654259
10 0.9457597 -5.6709834
11 -3.5905325 0.9457597
12 5.4411523 -3.5905325
13 -1.6672936 5.4411523
14 -2.7865613 -1.6672936
15 -2.2172448 -2.7865613
16 2.6559675 -2.2172448
17 3.0553994 2.6559675
18 3.2477355 3.0553994
19 0.8698169 3.2477355
20 2.4965039 0.8698169
21 -1.8119560 2.4965039
22 2.4236875 -1.8119560
23 -3.7143702 2.4236875
24 2.1756675 -3.7143702
25 -1.0449988 2.1756675
26 2.2965927 -1.0449988
27 1.0551433 2.2965927
28 -3.8594973 1.0551433
29 5.4767383 -3.8594973
30 3.2388077 5.4767383
31 -0.7632393 3.2388077
32 4.3024987 -0.7632393
33 6.8813701 4.3024987
34 -8.1602143 6.8813701
35 4.3601626 -8.1602143
36 0.1992655 4.3601626
37 0.7038718 0.1992655
38 -1.2786653 0.7038718
39 -0.6566598 -1.2786653
40 1.2206849 -0.6566598
41 -2.3473297 1.2206849
42 -7.2426091 -2.3473297
43 -4.1565198 -7.2426091
44 -4.7367314 -4.1565198
45 -0.5743614 -4.7367314
46 -3.9210982 -0.5743614
47 -1.5026892 -3.9210982
48 -3.7950658 -1.5026892
49 -1.5995584 -3.7950658
50 4.8372026 -1.5995584
51 0.1232347 4.8372026
52 -0.9261378 0.1232347
53 -2.3569004 -0.9261378
54 -2.2199110 -2.3569004
55 -0.8754832 -2.2199110
56 -3.6276971 -0.8754832
57 1.1759307 -3.6276971
58 8.7118653 1.1759307
59 4.4474293 8.7118653
60 -1.9574574 4.4474293
61 0.0750702 -1.9574574
62 4.0140281 0.0750702
63 1.3459570 4.0140281
64 NA 1.3459570
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 3.5329087 -2.0635621
[2,] -7.0825968 3.5329087
[3,] 0.3495696 -7.0825968
[4,] 0.9089827 0.3495696
[5,] -3.8279076 0.9089827
[6,] 2.9759768 -3.8279076
[7,] 4.9254253 2.9759768
[8,] 1.5654259 4.9254253
[9,] -5.6709834 1.5654259
[10,] 0.9457597 -5.6709834
[11,] -3.5905325 0.9457597
[12,] 5.4411523 -3.5905325
[13,] -1.6672936 5.4411523
[14,] -2.7865613 -1.6672936
[15,] -2.2172448 -2.7865613
[16,] 2.6559675 -2.2172448
[17,] 3.0553994 2.6559675
[18,] 3.2477355 3.0553994
[19,] 0.8698169 3.2477355
[20,] 2.4965039 0.8698169
[21,] -1.8119560 2.4965039
[22,] 2.4236875 -1.8119560
[23,] -3.7143702 2.4236875
[24,] 2.1756675 -3.7143702
[25,] -1.0449988 2.1756675
[26,] 2.2965927 -1.0449988
[27,] 1.0551433 2.2965927
[28,] -3.8594973 1.0551433
[29,] 5.4767383 -3.8594973
[30,] 3.2388077 5.4767383
[31,] -0.7632393 3.2388077
[32,] 4.3024987 -0.7632393
[33,] 6.8813701 4.3024987
[34,] -8.1602143 6.8813701
[35,] 4.3601626 -8.1602143
[36,] 0.1992655 4.3601626
[37,] 0.7038718 0.1992655
[38,] -1.2786653 0.7038718
[39,] -0.6566598 -1.2786653
[40,] 1.2206849 -0.6566598
[41,] -2.3473297 1.2206849
[42,] -7.2426091 -2.3473297
[43,] -4.1565198 -7.2426091
[44,] -4.7367314 -4.1565198
[45,] -0.5743614 -4.7367314
[46,] -3.9210982 -0.5743614
[47,] -1.5026892 -3.9210982
[48,] -3.7950658 -1.5026892
[49,] -1.5995584 -3.7950658
[50,] 4.8372026 -1.5995584
[51,] 0.1232347 4.8372026
[52,] -0.9261378 0.1232347
[53,] -2.3569004 -0.9261378
[54,] -2.2199110 -2.3569004
[55,] -0.8754832 -2.2199110
[56,] -3.6276971 -0.8754832
[57,] 1.1759307 -3.6276971
[58,] 8.7118653 1.1759307
[59,] 4.4474293 8.7118653
[60,] -1.9574574 4.4474293
[61,] 0.0750702 -1.9574574
[62,] 4.0140281 0.0750702
[63,] 1.3459570 4.0140281
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 3.5329087 -2.0635621
2 -7.0825968 3.5329087
3 0.3495696 -7.0825968
4 0.9089827 0.3495696
5 -3.8279076 0.9089827
6 2.9759768 -3.8279076
7 4.9254253 2.9759768
8 1.5654259 4.9254253
9 -5.6709834 1.5654259
10 0.9457597 -5.6709834
11 -3.5905325 0.9457597
12 5.4411523 -3.5905325
13 -1.6672936 5.4411523
14 -2.7865613 -1.6672936
15 -2.2172448 -2.7865613
16 2.6559675 -2.2172448
17 3.0553994 2.6559675
18 3.2477355 3.0553994
19 0.8698169 3.2477355
20 2.4965039 0.8698169
21 -1.8119560 2.4965039
22 2.4236875 -1.8119560
23 -3.7143702 2.4236875
24 2.1756675 -3.7143702
25 -1.0449988 2.1756675
26 2.2965927 -1.0449988
27 1.0551433 2.2965927
28 -3.8594973 1.0551433
29 5.4767383 -3.8594973
30 3.2388077 5.4767383
31 -0.7632393 3.2388077
32 4.3024987 -0.7632393
33 6.8813701 4.3024987
34 -8.1602143 6.8813701
35 4.3601626 -8.1602143
36 0.1992655 4.3601626
37 0.7038718 0.1992655
38 -1.2786653 0.7038718
39 -0.6566598 -1.2786653
40 1.2206849 -0.6566598
41 -2.3473297 1.2206849
42 -7.2426091 -2.3473297
43 -4.1565198 -7.2426091
44 -4.7367314 -4.1565198
45 -0.5743614 -4.7367314
46 -3.9210982 -0.5743614
47 -1.5026892 -3.9210982
48 -3.7950658 -1.5026892
49 -1.5995584 -3.7950658
50 4.8372026 -1.5995584
51 0.1232347 4.8372026
52 -0.9261378 0.1232347
53 -2.3569004 -0.9261378
54 -2.2199110 -2.3569004
55 -0.8754832 -2.2199110
56 -3.6276971 -0.8754832
57 1.1759307 -3.6276971
58 8.7118653 1.1759307
59 4.4474293 8.7118653
60 -1.9574574 4.4474293
61 0.0750702 -1.9574574
62 4.0140281 0.0750702
63 1.3459570 4.0140281
> 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/7mszm1291055503.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/8fjgp1291055503.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/9fjgp1291055503.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/10qafs1291055503.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/11bavx1291055503.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/12wtcl1291055503.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/13alsc1291055503.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/14w38i1291055503.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/15h4p61291055503.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/16dwnx1291055503.tab")
+ }
>
> try(system("convert tmp/1j90y1291055503.ps tmp/1j90y1291055503.png",intern=TRUE))
character(0)
> try(system("convert tmp/2j90y1291055503.ps tmp/2j90y1291055503.png",intern=TRUE))
character(0)
> try(system("convert tmp/3j90y1291055503.ps tmp/3j90y1291055503.png",intern=TRUE))
character(0)
> try(system("convert tmp/4t0hj1291055503.ps tmp/4t0hj1291055503.png",intern=TRUE))
character(0)
> try(system("convert tmp/5t0hj1291055503.ps tmp/5t0hj1291055503.png",intern=TRUE))
character(0)
> try(system("convert tmp/6mszm1291055503.ps tmp/6mszm1291055503.png",intern=TRUE))
character(0)
> try(system("convert tmp/7mszm1291055503.ps tmp/7mszm1291055503.png",intern=TRUE))
character(0)
> try(system("convert tmp/8fjgp1291055503.ps tmp/8fjgp1291055503.png",intern=TRUE))
character(0)
> try(system("convert tmp/9fjgp1291055503.ps tmp/9fjgp1291055503.png",intern=TRUE))
character(0)
> try(system("convert tmp/10qafs1291055503.ps tmp/10qafs1291055503.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.858 2.503 4.268