R version 2.9.0 (2009-04-17)
Copyright (C) 2009 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(25.4
+ ,23.4
+ ,11.5
+ ,9.9
+ ,8.1
+ ,27.9
+ ,25.4
+ ,23.4
+ ,11.5
+ ,9.9
+ ,26.1
+ ,27.9
+ ,25.4
+ ,23.4
+ ,11.5
+ ,18.8
+ ,26.1
+ ,27.9
+ ,25.4
+ ,23.4
+ ,14.1
+ ,18.8
+ ,26.1
+ ,27.9
+ ,25.4
+ ,11.5
+ ,14.1
+ ,18.8
+ ,26.1
+ ,27.9
+ ,15.8
+ ,11.5
+ ,14.1
+ ,18.8
+ ,26.1
+ ,12.4
+ ,15.8
+ ,11.5
+ ,14.1
+ ,18.8
+ ,4.5
+ ,12.4
+ ,15.8
+ ,11.5
+ ,14.1
+ ,-2.2
+ ,4.5
+ ,12.4
+ ,15.8
+ ,11.5
+ ,-4.2
+ ,-2.2
+ ,4.5
+ ,12.4
+ ,15.8
+ ,-9.4
+ ,-4.2
+ ,-2.2
+ ,4.5
+ ,12.4
+ ,-14.5
+ ,-9.4
+ ,-4.2
+ ,-2.2
+ ,4.5
+ ,-17.9
+ ,-14.5
+ ,-9.4
+ ,-4.2
+ ,-2.2
+ ,-15.1
+ ,-17.9
+ ,-14.5
+ ,-9.4
+ ,-4.2
+ ,-15.2
+ ,-15.1
+ ,-17.9
+ ,-14.5
+ ,-9.4
+ ,-15.7
+ ,-15.2
+ ,-15.1
+ ,-17.9
+ ,-14.5
+ ,-18
+ ,-15.7
+ ,-15.2
+ ,-15.1
+ ,-17.9
+ ,-18.1
+ ,-18
+ ,-15.7
+ ,-15.2
+ ,-15.1
+ ,-13.5
+ ,-18.1
+ ,-18
+ ,-15.7
+ ,-15.2
+ ,-9.9
+ ,-13.5
+ ,-18.1
+ ,-18
+ ,-15.7
+ ,-4.8
+ ,-9.9
+ ,-13.5
+ ,-18.1
+ ,-18
+ ,-1.7
+ ,-4.8
+ ,-9.9
+ ,-13.5
+ ,-18.1
+ ,-0.1
+ ,-1.7
+ ,-4.8
+ ,-9.9
+ ,-13.5
+ ,2.2
+ ,-0.1
+ ,-1.7
+ ,-4.8
+ ,-9.9
+ ,10.2
+ ,2.2
+ ,-0.1
+ ,-1.7
+ ,-4.8
+ ,7.6
+ ,10.2
+ ,2.2
+ ,-0.1
+ ,-1.7
+ ,10.8
+ ,7.6
+ ,10.2
+ ,2.2
+ ,-0.1
+ ,3.8
+ ,10.8
+ ,7.6
+ ,10.2
+ ,2.2
+ ,11
+ ,3.8
+ ,10.8
+ ,7.6
+ ,10.2
+ ,10.8
+ ,11
+ ,3.8
+ ,10.8
+ ,7.6
+ ,20.1
+ ,10.8
+ ,11
+ ,3.8
+ ,10.8
+ ,14.9
+ ,20.1
+ ,10.8
+ ,11
+ ,3.8
+ ,13
+ ,14.9
+ ,20.1
+ ,10.8
+ ,11
+ ,10.9
+ ,13
+ ,14.9
+ ,20.1
+ ,10.8
+ ,9.6
+ ,10.9
+ ,13
+ ,14.9
+ ,20.1
+ ,4
+ ,9.6
+ ,10.9
+ ,13
+ ,14.9
+ ,-1.1
+ ,4
+ ,9.6
+ ,10.9
+ ,13
+ ,-7.7
+ ,-1.1
+ ,4
+ ,9.6
+ ,10.9
+ ,-8.9
+ ,-7.7
+ ,-1.1
+ ,4
+ ,9.6
+ ,-8
+ ,-8.9
+ ,-7.7
+ ,-1.1
+ ,4
+ ,-7.1
+ ,-8
+ ,-8.9
+ ,-7.7
+ ,-1.1
+ ,-5.3
+ ,-7.1
+ ,-8
+ ,-8.9
+ ,-7.7
+ ,-2.5
+ ,-5.3
+ ,-7.1
+ ,-8
+ ,-8.9
+ ,-2.4
+ ,-2.5
+ ,-5.3
+ ,-7.1
+ ,-8
+ ,-2.9
+ ,-2.4
+ ,-2.5
+ ,-5.3
+ ,-7.1
+ ,-4.8
+ ,-2.9
+ ,-2.4
+ ,-2.5
+ ,-5.3
+ ,-7.2
+ ,-4.8
+ ,-2.9
+ ,-2.4
+ ,-2.5
+ ,1.7
+ ,-7.2
+ ,-4.8
+ ,-2.9
+ ,-2.4
+ ,2.2
+ ,1.7
+ ,-7.2
+ ,-4.8
+ ,-2.9
+ ,13.4
+ ,2.2
+ ,1.7
+ ,-7.2
+ ,-4.8
+ ,12.3
+ ,13.4
+ ,2.2
+ ,1.7
+ ,-7.2
+ ,13.7
+ ,12.3
+ ,13.4
+ ,2.2
+ ,1.7
+ ,4.4
+ ,13.7
+ ,12.3
+ ,13.4
+ ,2.2
+ ,-2.5
+ ,4.4
+ ,13.7
+ ,12.3
+ ,13.4)
+ ,dim=c(5
+ ,55)
+ ,dimnames=list(c('Ye'
+ ,'Ye-1'
+ ,'Ye-2'
+ ,'Ye-3'
+ ,'Ye-4')
+ ,1:55))
> y <- array(NA,dim=c(5,55),dimnames=list(c('Ye','Ye-1','Ye-2','Ye-3','Ye-4'),1:55))
> 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
Ye Ye-1 Ye-2 Ye-3 Ye-4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 25.4 23.4 11.5 9.9 8.1 1 0 0 0 0 0 0 0 0 0 0 1
2 27.9 25.4 23.4 11.5 9.9 0 1 0 0 0 0 0 0 0 0 0 2
3 26.1 27.9 25.4 23.4 11.5 0 0 1 0 0 0 0 0 0 0 0 3
4 18.8 26.1 27.9 25.4 23.4 0 0 0 1 0 0 0 0 0 0 0 4
5 14.1 18.8 26.1 27.9 25.4 0 0 0 0 1 0 0 0 0 0 0 5
6 11.5 14.1 18.8 26.1 27.9 0 0 0 0 0 1 0 0 0 0 0 6
7 15.8 11.5 14.1 18.8 26.1 0 0 0 0 0 0 1 0 0 0 0 7
8 12.4 15.8 11.5 14.1 18.8 0 0 0 0 0 0 0 1 0 0 0 8
9 4.5 12.4 15.8 11.5 14.1 0 0 0 0 0 0 0 0 1 0 0 9
10 -2.2 4.5 12.4 15.8 11.5 0 0 0 0 0 0 0 0 0 1 0 10
11 -4.2 -2.2 4.5 12.4 15.8 0 0 0 0 0 0 0 0 0 0 1 11
12 -9.4 -4.2 -2.2 4.5 12.4 0 0 0 0 0 0 0 0 0 0 0 12
13 -14.5 -9.4 -4.2 -2.2 4.5 1 0 0 0 0 0 0 0 0 0 0 13
14 -17.9 -14.5 -9.4 -4.2 -2.2 0 1 0 0 0 0 0 0 0 0 0 14
15 -15.1 -17.9 -14.5 -9.4 -4.2 0 0 1 0 0 0 0 0 0 0 0 15
16 -15.2 -15.1 -17.9 -14.5 -9.4 0 0 0 1 0 0 0 0 0 0 0 16
17 -15.7 -15.2 -15.1 -17.9 -14.5 0 0 0 0 1 0 0 0 0 0 0 17
18 -18.0 -15.7 -15.2 -15.1 -17.9 0 0 0 0 0 1 0 0 0 0 0 18
19 -18.1 -18.0 -15.7 -15.2 -15.1 0 0 0 0 0 0 1 0 0 0 0 19
20 -13.5 -18.1 -18.0 -15.7 -15.2 0 0 0 0 0 0 0 1 0 0 0 20
21 -9.9 -13.5 -18.1 -18.0 -15.7 0 0 0 0 0 0 0 0 1 0 0 21
22 -4.8 -9.9 -13.5 -18.1 -18.0 0 0 0 0 0 0 0 0 0 1 0 22
23 -1.7 -4.8 -9.9 -13.5 -18.1 0 0 0 0 0 0 0 0 0 0 1 23
24 -0.1 -1.7 -4.8 -9.9 -13.5 0 0 0 0 0 0 0 0 0 0 0 24
25 2.2 -0.1 -1.7 -4.8 -9.9 1 0 0 0 0 0 0 0 0 0 0 25
26 10.2 2.2 -0.1 -1.7 -4.8 0 1 0 0 0 0 0 0 0 0 0 26
27 7.6 10.2 2.2 -0.1 -1.7 0 0 1 0 0 0 0 0 0 0 0 27
28 10.8 7.6 10.2 2.2 -0.1 0 0 0 1 0 0 0 0 0 0 0 28
29 3.8 10.8 7.6 10.2 2.2 0 0 0 0 1 0 0 0 0 0 0 29
30 11.0 3.8 10.8 7.6 10.2 0 0 0 0 0 1 0 0 0 0 0 30
31 10.8 11.0 3.8 10.8 7.6 0 0 0 0 0 0 1 0 0 0 0 31
32 20.1 10.8 11.0 3.8 10.8 0 0 0 0 0 0 0 1 0 0 0 32
33 14.9 20.1 10.8 11.0 3.8 0 0 0 0 0 0 0 0 1 0 0 33
34 13.0 14.9 20.1 10.8 11.0 0 0 0 0 0 0 0 0 0 1 0 34
35 10.9 13.0 14.9 20.1 10.8 0 0 0 0 0 0 0 0 0 0 1 35
36 9.6 10.9 13.0 14.9 20.1 0 0 0 0 0 0 0 0 0 0 0 36
37 4.0 9.6 10.9 13.0 14.9 1 0 0 0 0 0 0 0 0 0 0 37
38 -1.1 4.0 9.6 10.9 13.0 0 1 0 0 0 0 0 0 0 0 0 38
39 -7.7 -1.1 4.0 9.6 10.9 0 0 1 0 0 0 0 0 0 0 0 39
40 -8.9 -7.7 -1.1 4.0 9.6 0 0 0 1 0 0 0 0 0 0 0 40
41 -8.0 -8.9 -7.7 -1.1 4.0 0 0 0 0 1 0 0 0 0 0 0 41
42 -7.1 -8.0 -8.9 -7.7 -1.1 0 0 0 0 0 1 0 0 0 0 0 42
43 -5.3 -7.1 -8.0 -8.9 -7.7 0 0 0 0 0 0 1 0 0 0 0 43
44 -2.5 -5.3 -7.1 -8.0 -8.9 0 0 0 0 0 0 0 1 0 0 0 44
45 -2.4 -2.5 -5.3 -7.1 -8.0 0 0 0 0 0 0 0 0 1 0 0 45
46 -2.9 -2.4 -2.5 -5.3 -7.1 0 0 0 0 0 0 0 0 0 1 0 46
47 -4.8 -2.9 -2.4 -2.5 -5.3 0 0 0 0 0 0 0 0 0 0 1 47
48 -7.2 -4.8 -2.9 -2.4 -2.5 0 0 0 0 0 0 0 0 0 0 0 48
49 1.7 -7.2 -4.8 -2.9 -2.4 1 0 0 0 0 0 0 0 0 0 0 49
50 2.2 1.7 -7.2 -4.8 -2.9 0 1 0 0 0 0 0 0 0 0 0 50
51 13.4 2.2 1.7 -7.2 -4.8 0 0 1 0 0 0 0 0 0 0 0 51
52 12.3 13.4 2.2 1.7 -7.2 0 0 0 1 0 0 0 0 0 0 0 52
53 13.7 12.3 13.4 2.2 1.7 0 0 0 0 1 0 0 0 0 0 0 53
54 4.4 13.7 12.3 13.4 2.2 0 0 0 0 0 1 0 0 0 0 0 54
55 -2.5 4.4 13.7 12.3 13.4 0 0 0 0 0 0 1 0 0 0 0 55
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) `Ye-1` `Ye-2` `Ye-3` `Ye-4` M1
-1.923021 1.170772 0.192986 -0.758568 0.279151 2.425975
M2 M3 M4 M5 M6 M7
2.084145 2.768231 0.787245 0.448670 2.108750 2.541080
M8 M9 M10 M11 t
3.636192 -1.495497 0.486247 3.498512 0.004496
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-6.1706 -2.1259 -0.4162 2.1218 8.8027
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.923021 2.384285 -0.807 0.424951
`Ye-1` 1.170772 0.145955 8.021 1.07e-09 ***
`Ye-2` 0.192986 0.197242 0.978 0.334052
`Ye-3` -0.758568 0.207656 -3.653 0.000779 ***
`Ye-4` 0.279151 0.154243 1.810 0.078241 .
M1 2.425975 2.758840 0.879 0.384741
M2 2.084145 2.764157 0.754 0.455503
M3 2.768231 2.776054 0.997 0.324984
M4 0.787245 2.766462 0.285 0.777521
M5 0.448670 2.777923 0.162 0.872545
M6 2.108750 2.787292 0.757 0.453979
M7 2.541080 2.748009 0.925 0.360958
M8 3.636192 2.918006 1.246 0.220347
M9 -1.495497 2.948789 -0.507 0.614976
M10 0.486247 3.019172 0.161 0.872905
M11 3.498512 3.024528 1.157 0.254610
t 0.004496 0.036684 0.123 0.903100
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.081 on 38 degrees of freedom
Multiple R-squared: 0.9186, Adjusted R-squared: 0.8843
F-statistic: 26.8 on 16 and 38 DF, p-value: 7.503e-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.3603609 0.7207218 0.6396391
[2,] 0.6314247 0.7371506 0.3685753
[3,] 0.4749426 0.9498853 0.5250574
[4,] 0.3663544 0.7327088 0.6336456
[5,] 0.2805248 0.5610497 0.7194752
[6,] 0.1986895 0.3973790 0.8013105
[7,] 0.1501046 0.3002092 0.8498954
[8,] 0.4357718 0.8715436 0.5642282
[9,] 0.3976656 0.7953312 0.6023344
[10,] 0.4362288 0.8724576 0.5637712
[11,] 0.4671470 0.9342940 0.5328530
[12,] 0.3909965 0.7819930 0.6090035
[13,] 0.3138385 0.6276769 0.6861615
[14,] 0.2044973 0.4089945 0.7955027
[15,] 0.1364913 0.2729826 0.8635087
[16,] 0.2245318 0.4490636 0.7754682
> postscript(file="/var/www/html/rcomp/tmp/1o0jy1292958655.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/rcomp/tmp/2yr111292958655.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/rcomp/tmp/3yr111292958655.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/rcomp/tmp/4yr111292958655.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/rcomp/tmp/5r00m1292958655.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 = 55
Frequency = 1
1 2 3 4 5 6 7
0.5258460 -0.5636604 2.2151762 -3.2881644 2.5780462 3.1615981 5.9407375
8 9 10 11 12 13 14
-4.6188985 -4.9011908 0.3054585 0.8779817 -2.2370290 -6.1706305 -1.9056569
15 16 17 18 19 20 21
1.7843609 -1.3782728 -3.1229419 -3.4097302 -2.0147652 1.6952031 3.4510099
22 23 24 25 26 27 28
2.0284441 -1.0366756 0.8904698 1.1522622 7.4159367 -5.3343495 2.6403261
29 30 31 32 33 34 35
-1.8438018 7.0639886 2.5017166 3.3435025 -0.1631430 -1.9177362 3.3040136
36 37 38 39 40 41 42
1.7826681 -4.3102251 -3.3282913 -3.9651361 1.6375936 3.2448528 -1.9247182
43 44 45 46 47 48 49
-0.8568139 -0.4198071 1.6133239 -0.4161664 -3.1453197 -0.4361089 8.8027474
50 51 52 53 54 55
-1.6183281 5.2999486 0.3885175 -0.8561554 -4.8911382 -5.5708752
> postscript(file="/var/www/html/rcomp/tmp/6r00m1292958655.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 = 55
Frequency = 1
lag(myerror, k = 1) myerror
0 0.5258460 NA
1 -0.5636604 0.5258460
2 2.2151762 -0.5636604
3 -3.2881644 2.2151762
4 2.5780462 -3.2881644
5 3.1615981 2.5780462
6 5.9407375 3.1615981
7 -4.6188985 5.9407375
8 -4.9011908 -4.6188985
9 0.3054585 -4.9011908
10 0.8779817 0.3054585
11 -2.2370290 0.8779817
12 -6.1706305 -2.2370290
13 -1.9056569 -6.1706305
14 1.7843609 -1.9056569
15 -1.3782728 1.7843609
16 -3.1229419 -1.3782728
17 -3.4097302 -3.1229419
18 -2.0147652 -3.4097302
19 1.6952031 -2.0147652
20 3.4510099 1.6952031
21 2.0284441 3.4510099
22 -1.0366756 2.0284441
23 0.8904698 -1.0366756
24 1.1522622 0.8904698
25 7.4159367 1.1522622
26 -5.3343495 7.4159367
27 2.6403261 -5.3343495
28 -1.8438018 2.6403261
29 7.0639886 -1.8438018
30 2.5017166 7.0639886
31 3.3435025 2.5017166
32 -0.1631430 3.3435025
33 -1.9177362 -0.1631430
34 3.3040136 -1.9177362
35 1.7826681 3.3040136
36 -4.3102251 1.7826681
37 -3.3282913 -4.3102251
38 -3.9651361 -3.3282913
39 1.6375936 -3.9651361
40 3.2448528 1.6375936
41 -1.9247182 3.2448528
42 -0.8568139 -1.9247182
43 -0.4198071 -0.8568139
44 1.6133239 -0.4198071
45 -0.4161664 1.6133239
46 -3.1453197 -0.4161664
47 -0.4361089 -3.1453197
48 8.8027474 -0.4361089
49 -1.6183281 8.8027474
50 5.2999486 -1.6183281
51 0.3885175 5.2999486
52 -0.8561554 0.3885175
53 -4.8911382 -0.8561554
54 -5.5708752 -4.8911382
55 NA -5.5708752
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.5636604 0.5258460
[2,] 2.2151762 -0.5636604
[3,] -3.2881644 2.2151762
[4,] 2.5780462 -3.2881644
[5,] 3.1615981 2.5780462
[6,] 5.9407375 3.1615981
[7,] -4.6188985 5.9407375
[8,] -4.9011908 -4.6188985
[9,] 0.3054585 -4.9011908
[10,] 0.8779817 0.3054585
[11,] -2.2370290 0.8779817
[12,] -6.1706305 -2.2370290
[13,] -1.9056569 -6.1706305
[14,] 1.7843609 -1.9056569
[15,] -1.3782728 1.7843609
[16,] -3.1229419 -1.3782728
[17,] -3.4097302 -3.1229419
[18,] -2.0147652 -3.4097302
[19,] 1.6952031 -2.0147652
[20,] 3.4510099 1.6952031
[21,] 2.0284441 3.4510099
[22,] -1.0366756 2.0284441
[23,] 0.8904698 -1.0366756
[24,] 1.1522622 0.8904698
[25,] 7.4159367 1.1522622
[26,] -5.3343495 7.4159367
[27,] 2.6403261 -5.3343495
[28,] -1.8438018 2.6403261
[29,] 7.0639886 -1.8438018
[30,] 2.5017166 7.0639886
[31,] 3.3435025 2.5017166
[32,] -0.1631430 3.3435025
[33,] -1.9177362 -0.1631430
[34,] 3.3040136 -1.9177362
[35,] 1.7826681 3.3040136
[36,] -4.3102251 1.7826681
[37,] -3.3282913 -4.3102251
[38,] -3.9651361 -3.3282913
[39,] 1.6375936 -3.9651361
[40,] 3.2448528 1.6375936
[41,] -1.9247182 3.2448528
[42,] -0.8568139 -1.9247182
[43,] -0.4198071 -0.8568139
[44,] 1.6133239 -0.4198071
[45,] -0.4161664 1.6133239
[46,] -3.1453197 -0.4161664
[47,] -0.4361089 -3.1453197
[48,] 8.8027474 -0.4361089
[49,] -1.6183281 8.8027474
[50,] 5.2999486 -1.6183281
[51,] 0.3885175 5.2999486
[52,] -0.8561554 0.3885175
[53,] -4.8911382 -0.8561554
[54,] -5.5708752 -4.8911382
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.5636604 0.5258460
2 2.2151762 -0.5636604
3 -3.2881644 2.2151762
4 2.5780462 -3.2881644
5 3.1615981 2.5780462
6 5.9407375 3.1615981
7 -4.6188985 5.9407375
8 -4.9011908 -4.6188985
9 0.3054585 -4.9011908
10 0.8779817 0.3054585
11 -2.2370290 0.8779817
12 -6.1706305 -2.2370290
13 -1.9056569 -6.1706305
14 1.7843609 -1.9056569
15 -1.3782728 1.7843609
16 -3.1229419 -1.3782728
17 -3.4097302 -3.1229419
18 -2.0147652 -3.4097302
19 1.6952031 -2.0147652
20 3.4510099 1.6952031
21 2.0284441 3.4510099
22 -1.0366756 2.0284441
23 0.8904698 -1.0366756
24 1.1522622 0.8904698
25 7.4159367 1.1522622
26 -5.3343495 7.4159367
27 2.6403261 -5.3343495
28 -1.8438018 2.6403261
29 7.0639886 -1.8438018
30 2.5017166 7.0639886
31 3.3435025 2.5017166
32 -0.1631430 3.3435025
33 -1.9177362 -0.1631430
34 3.3040136 -1.9177362
35 1.7826681 3.3040136
36 -4.3102251 1.7826681
37 -3.3282913 -4.3102251
38 -3.9651361 -3.3282913
39 1.6375936 -3.9651361
40 3.2448528 1.6375936
41 -1.9247182 3.2448528
42 -0.8568139 -1.9247182
43 -0.4198071 -0.8568139
44 1.6133239 -0.4198071
45 -0.4161664 1.6133239
46 -3.1453197 -0.4161664
47 -0.4361089 -3.1453197
48 8.8027474 -0.4361089
49 -1.6183281 8.8027474
50 5.2999486 -1.6183281
51 0.3885175 5.2999486
52 -0.8561554 0.3885175
53 -4.8911382 -0.8561554
54 -5.5708752 -4.8911382
> 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/rcomp/tmp/72shp1292958655.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/rcomp/tmp/82shp1292958655.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/rcomp/tmp/9cjys1292958655.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/rcomp/tmp/10nsyv1292958655.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/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/rcomp/tmp/11yjff1292958655.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/rcomp/tmp/121kd31292958655.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/rcomp/tmp/138lsf1292958655.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/rcomp/tmp/14b3r31292958655.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/rcomp/tmp/15fm8r1292958655.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/rcomp/tmp/160m6f1292958655.tab")
+ }
> try(system("convert tmp/1o0jy1292958655.ps tmp/1o0jy1292958655.png",intern=TRUE))
character(0)
> try(system("convert tmp/2yr111292958655.ps tmp/2yr111292958655.png",intern=TRUE))
character(0)
> try(system("convert tmp/3yr111292958655.ps tmp/3yr111292958655.png",intern=TRUE))
character(0)
> try(system("convert tmp/4yr111292958655.ps tmp/4yr111292958655.png",intern=TRUE))
character(0)
> try(system("convert tmp/5r00m1292958655.ps tmp/5r00m1292958655.png",intern=TRUE))
character(0)
> try(system("convert tmp/6r00m1292958655.ps tmp/6r00m1292958655.png",intern=TRUE))
character(0)
> try(system("convert tmp/72shp1292958655.ps tmp/72shp1292958655.png",intern=TRUE))
character(0)
> try(system("convert tmp/82shp1292958655.ps tmp/82shp1292958655.png",intern=TRUE))
character(0)
> try(system("convert tmp/9cjys1292958655.ps tmp/9cjys1292958655.png",intern=TRUE))
character(0)
> try(system("convert tmp/10nsyv1292958655.ps tmp/10nsyv1292958655.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.361 1.639 7.051