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(112.3
+ ,1
+ ,117.2
+ ,96.8
+ ,80
+ ,126.1
+ ,117.3
+ ,1
+ ,112.3
+ ,117.2
+ ,96.8
+ ,80
+ ,111.1
+ ,1
+ ,117.3
+ ,112.3
+ ,117.2
+ ,96.8
+ ,102.2
+ ,1
+ ,111.1
+ ,117.3
+ ,112.3
+ ,117.2
+ ,104.3
+ ,1
+ ,102.2
+ ,111.1
+ ,117.3
+ ,112.3
+ ,122.9
+ ,0
+ ,104.3
+ ,102.2
+ ,111.1
+ ,117.3
+ ,107.6
+ ,0
+ ,122.9
+ ,104.3
+ ,102.2
+ ,111.1
+ ,121.3
+ ,0
+ ,107.6
+ ,122.9
+ ,104.3
+ ,102.2
+ ,131.5
+ ,0
+ ,121.3
+ ,107.6
+ ,122.9
+ ,104.3
+ ,89
+ ,0
+ ,131.5
+ ,121.3
+ ,107.6
+ ,122.9
+ ,104.4
+ ,0
+ ,89
+ ,131.5
+ ,121.3
+ ,107.6
+ ,128.9
+ ,0
+ ,104.4
+ ,89
+ ,131.5
+ ,121.3
+ ,135.9
+ ,0
+ ,128.9
+ ,104.4
+ ,89
+ ,131.5
+ ,133.3
+ ,0
+ ,135.9
+ ,128.9
+ ,104.4
+ ,89
+ ,121.3
+ ,0
+ ,133.3
+ ,135.9
+ ,128.9
+ ,104.4
+ ,120.5
+ ,0
+ ,121.3
+ ,133.3
+ ,135.9
+ ,128.9
+ ,120.4
+ ,0
+ ,120.5
+ ,121.3
+ ,133.3
+ ,135.9
+ ,137.9
+ ,0
+ ,120.4
+ ,120.5
+ ,121.3
+ ,133.3
+ ,126.1
+ ,0
+ ,137.9
+ ,120.4
+ ,120.5
+ ,121.3
+ ,133.2
+ ,0
+ ,126.1
+ ,137.9
+ ,120.4
+ ,120.5
+ ,151.1
+ ,0
+ ,133.2
+ ,126.1
+ ,137.9
+ ,120.4
+ ,105
+ ,0
+ ,151.1
+ ,133.2
+ ,126.1
+ ,137.9
+ ,119
+ ,0
+ ,105
+ ,151.1
+ ,133.2
+ ,126.1
+ ,140.4
+ ,0
+ ,119
+ ,105
+ ,151.1
+ ,133.2
+ ,156.6
+ ,1
+ ,140.4
+ ,119
+ ,105
+ ,151.1
+ ,137.1
+ ,1
+ ,156.6
+ ,140.4
+ ,119
+ ,105
+ ,122.7
+ ,1
+ ,137.1
+ ,156.6
+ ,140.4
+ ,119
+ ,125.8
+ ,1
+ ,122.7
+ ,137.1
+ ,156.6
+ ,140.4
+ ,139.3
+ ,1
+ ,125.8
+ ,122.7
+ ,137.1
+ ,156.6
+ ,134.9
+ ,1
+ ,139.3
+ ,125.8
+ ,122.7
+ ,137.1
+ ,149.2
+ ,1
+ ,134.9
+ ,139.3
+ ,125.8
+ ,122.7
+ ,132.3
+ ,1
+ ,149.2
+ ,134.9
+ ,139.3
+ ,125.8
+ ,149
+ ,1
+ ,132.3
+ ,149.2
+ ,134.9
+ ,139.3
+ ,117.2
+ ,1
+ ,149
+ ,132.3
+ ,149.2
+ ,134.9
+ ,119.6
+ ,1
+ ,117.2
+ ,149
+ ,132.3
+ ,149.2
+ ,152
+ ,1
+ ,119.6
+ ,117.2
+ ,149
+ ,132.3
+ ,149.4
+ ,1
+ ,152
+ ,119.6
+ ,117.2
+ ,149
+ ,127.3
+ ,1
+ ,149.4
+ ,152
+ ,119.6
+ ,117.2
+ ,114.1
+ ,1
+ ,127.3
+ ,149.4
+ ,152
+ ,119.6
+ ,102.1
+ ,1
+ ,114.1
+ ,127.3
+ ,149.4
+ ,152
+ ,107.7
+ ,1
+ ,102.1
+ ,114.1
+ ,127.3
+ ,149.4
+ ,104.4
+ ,1
+ ,107.7
+ ,102.1
+ ,114.1
+ ,127.3
+ ,102.1
+ ,1
+ ,104.4
+ ,107.7
+ ,102.1
+ ,114.1
+ ,96
+ ,1
+ ,102.1
+ ,104.4
+ ,107.7
+ ,102.1
+ ,109.3
+ ,1
+ ,96
+ ,102.1
+ ,104.4
+ ,107.7
+ ,90
+ ,1
+ ,109.3
+ ,96
+ ,102.1
+ ,104.4
+ ,83.9
+ ,1
+ ,90
+ ,109.3
+ ,96
+ ,102.1
+ ,112
+ ,1
+ ,83.9
+ ,90
+ ,109.3
+ ,96
+ ,114.3
+ ,1
+ ,112
+ ,83.9
+ ,90
+ ,109.3
+ ,103.6
+ ,1
+ ,114.3
+ ,112
+ ,83.9
+ ,90
+ ,91.7
+ ,1
+ ,103.6
+ ,114.3
+ ,112
+ ,83.9
+ ,80.8
+ ,1
+ ,91.7
+ ,103.6
+ ,114.3
+ ,112
+ ,87.2
+ ,1
+ ,80.8
+ ,91.7
+ ,103.6
+ ,114.3
+ ,109.2
+ ,1
+ ,87.2
+ ,80.8
+ ,91.7
+ ,103.6
+ ,102.7
+ ,1
+ ,109.2
+ ,87.2
+ ,80.8
+ ,91.7
+ ,95.1
+ ,1
+ ,102.7
+ ,109.2
+ ,87.2
+ ,80.8
+ ,117.5
+ ,1
+ ,95.1
+ ,102.7
+ ,109.2
+ ,87.2
+ ,85.1
+ ,1
+ ,117.5
+ ,95.1
+ ,102.7
+ ,109.2
+ ,92.1
+ ,1
+ ,85.1
+ ,117.5
+ ,95.1
+ ,102.7
+ ,113.5
+ ,1
+ ,92.1
+ ,85.1
+ ,117.5
+ ,95.1)
+ ,dim=c(6
+ ,60)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y1'
+ ,'Y2'
+ ,'Y3'
+ ,'Y4')
+ ,1:60))
> y <- array(NA,dim=c(6,60),dimnames=list(c('Y','X','Y1','Y2','Y3','Y4'),1:60))
> 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 Y1 Y2 Y3 Y4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 112.3 1 117.2 96.8 80.0 126.1 1 0 0 0 0 0 0 0 0 0 0 1
2 117.3 1 112.3 117.2 96.8 80.0 0 1 0 0 0 0 0 0 0 0 0 2
3 111.1 1 117.3 112.3 117.2 96.8 0 0 1 0 0 0 0 0 0 0 0 3
4 102.2 1 111.1 117.3 112.3 117.2 0 0 0 1 0 0 0 0 0 0 0 4
5 104.3 1 102.2 111.1 117.3 112.3 0 0 0 0 1 0 0 0 0 0 0 5
6 122.9 0 104.3 102.2 111.1 117.3 0 0 0 0 0 1 0 0 0 0 0 6
7 107.6 0 122.9 104.3 102.2 111.1 0 0 0 0 0 0 1 0 0 0 0 7
8 121.3 0 107.6 122.9 104.3 102.2 0 0 0 0 0 0 0 1 0 0 0 8
9 131.5 0 121.3 107.6 122.9 104.3 0 0 0 0 0 0 0 0 1 0 0 9
10 89.0 0 131.5 121.3 107.6 122.9 0 0 0 0 0 0 0 0 0 1 0 10
11 104.4 0 89.0 131.5 121.3 107.6 0 0 0 0 0 0 0 0 0 0 1 11
12 128.9 0 104.4 89.0 131.5 121.3 0 0 0 0 0 0 0 0 0 0 0 12
13 135.9 0 128.9 104.4 89.0 131.5 1 0 0 0 0 0 0 0 0 0 0 13
14 133.3 0 135.9 128.9 104.4 89.0 0 1 0 0 0 0 0 0 0 0 0 14
15 121.3 0 133.3 135.9 128.9 104.4 0 0 1 0 0 0 0 0 0 0 0 15
16 120.5 0 121.3 133.3 135.9 128.9 0 0 0 1 0 0 0 0 0 0 0 16
17 120.4 0 120.5 121.3 133.3 135.9 0 0 0 0 1 0 0 0 0 0 0 17
18 137.9 0 120.4 120.5 121.3 133.3 0 0 0 0 0 1 0 0 0 0 0 18
19 126.1 0 137.9 120.4 120.5 121.3 0 0 0 0 0 0 1 0 0 0 0 19
20 133.2 0 126.1 137.9 120.4 120.5 0 0 0 0 0 0 0 1 0 0 0 20
21 151.1 0 133.2 126.1 137.9 120.4 0 0 0 0 0 0 0 0 1 0 0 21
22 105.0 0 151.1 133.2 126.1 137.9 0 0 0 0 0 0 0 0 0 1 0 22
23 119.0 0 105.0 151.1 133.2 126.1 0 0 0 0 0 0 0 0 0 0 1 23
24 140.4 0 119.0 105.0 151.1 133.2 0 0 0 0 0 0 0 0 0 0 0 24
25 156.6 1 140.4 119.0 105.0 151.1 1 0 0 0 0 0 0 0 0 0 0 25
26 137.1 1 156.6 140.4 119.0 105.0 0 1 0 0 0 0 0 0 0 0 0 26
27 122.7 1 137.1 156.6 140.4 119.0 0 0 1 0 0 0 0 0 0 0 0 27
28 125.8 1 122.7 137.1 156.6 140.4 0 0 0 1 0 0 0 0 0 0 0 28
29 139.3 1 125.8 122.7 137.1 156.6 0 0 0 0 1 0 0 0 0 0 0 29
30 134.9 1 139.3 125.8 122.7 137.1 0 0 0 0 0 1 0 0 0 0 0 30
31 149.2 1 134.9 139.3 125.8 122.7 0 0 0 0 0 0 1 0 0 0 0 31
32 132.3 1 149.2 134.9 139.3 125.8 0 0 0 0 0 0 0 1 0 0 0 32
33 149.0 1 132.3 149.2 134.9 139.3 0 0 0 0 0 0 0 0 1 0 0 33
34 117.2 1 149.0 132.3 149.2 134.9 0 0 0 0 0 0 0 0 0 1 0 34
35 119.6 1 117.2 149.0 132.3 149.2 0 0 0 0 0 0 0 0 0 0 1 35
36 152.0 1 119.6 117.2 149.0 132.3 0 0 0 0 0 0 0 0 0 0 0 36
37 149.4 1 152.0 119.6 117.2 149.0 1 0 0 0 0 0 0 0 0 0 0 37
38 127.3 1 149.4 152.0 119.6 117.2 0 1 0 0 0 0 0 0 0 0 0 38
39 114.1 1 127.3 149.4 152.0 119.6 0 0 1 0 0 0 0 0 0 0 0 39
40 102.1 1 114.1 127.3 149.4 152.0 0 0 0 1 0 0 0 0 0 0 0 40
41 107.7 1 102.1 114.1 127.3 149.4 0 0 0 0 1 0 0 0 0 0 0 41
42 104.4 1 107.7 102.1 114.1 127.3 0 0 0 0 0 1 0 0 0 0 0 42
43 102.1 1 104.4 107.7 102.1 114.1 0 0 0 0 0 0 1 0 0 0 0 43
44 96.0 1 102.1 104.4 107.7 102.1 0 0 0 0 0 0 0 1 0 0 0 44
45 109.3 1 96.0 102.1 104.4 107.7 0 0 0 0 0 0 0 0 1 0 0 45
46 90.0 1 109.3 96.0 102.1 104.4 0 0 0 0 0 0 0 0 0 1 0 46
47 83.9 1 90.0 109.3 96.0 102.1 0 0 0 0 0 0 0 0 0 0 1 47
48 112.0 1 83.9 90.0 109.3 96.0 0 0 0 0 0 0 0 0 0 0 0 48
49 114.3 1 112.0 83.9 90.0 109.3 1 0 0 0 0 0 0 0 0 0 0 49
50 103.6 1 114.3 112.0 83.9 90.0 0 1 0 0 0 0 0 0 0 0 0 50
51 91.7 1 103.6 114.3 112.0 83.9 0 0 1 0 0 0 0 0 0 0 0 51
52 80.8 1 91.7 103.6 114.3 112.0 0 0 0 1 0 0 0 0 0 0 0 52
53 87.2 1 80.8 91.7 103.6 114.3 0 0 0 0 1 0 0 0 0 0 0 53
54 109.2 1 87.2 80.8 91.7 103.6 0 0 0 0 0 1 0 0 0 0 0 54
55 102.7 1 109.2 87.2 80.8 91.7 0 0 0 0 0 0 1 0 0 0 0 55
56 95.1 1 102.7 109.2 87.2 80.8 0 0 0 0 0 0 0 1 0 0 0 56
57 117.5 1 95.1 102.7 109.2 87.2 0 0 0 0 0 0 0 0 1 0 0 57
58 85.1 1 117.5 95.1 102.7 109.2 0 0 0 0 0 0 0 0 0 1 0 58
59 92.1 1 85.1 117.5 95.1 102.7 0 0 0 0 0 0 0 0 0 0 1 59
60 113.5 1 92.1 85.1 117.5 95.1 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X Y1 Y2 Y3 Y4
36.8237 -0.8244 0.3761 0.3528 0.2843 -0.1245
M1 M2 M3 M4 M5 M6
3.2357 -23.9939 -39.1094 -34.9214 -19.7511 -7.6898
M7 M8 M9 M10 M11 t
-17.4063 -23.5023 -7.2285 -44.4355 -30.5172 -0.0939
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-13.0944 -4.2743 -0.7317 4.6872 15.6773
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 36.82370 10.66629 3.452 0.001281 **
X -0.82445 3.01600 -0.273 0.785918
Y1 0.37610 0.15198 2.475 0.017460 *
Y2 0.35284 0.15702 2.247 0.029950 *
Y3 0.28433 0.15340 1.854 0.070841 .
Y4 -0.12446 0.14806 -0.841 0.405347
M1 3.23567 9.78102 0.331 0.742432
M2 -23.99390 10.61758 -2.260 0.029079 *
M3 -39.10940 8.15636 -4.795 2.07e-05 ***
M4 -34.92136 5.99615 -5.824 7.11e-07 ***
M5 -19.75110 6.16858 -3.202 0.002603 **
M6 -7.68983 6.50772 -1.182 0.243993
M7 -17.40631 7.73829 -2.249 0.029790 *
M8 -23.50234 7.86672 -2.988 0.004680 **
M9 -7.22853 6.46008 -1.119 0.269519
M10 -44.43546 7.41513 -5.993 4.07e-07 ***
M11 -30.51724 8.37822 -3.642 0.000736 ***
t -0.09391 0.08445 -1.112 0.272479
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 7.791 on 42 degrees of freedom
Multiple R-squared: 0.8799, Adjusted R-squared: 0.8313
F-statistic: 18.1 on 17 and 42 DF, p-value: 4.269e-14
> 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.0220619366 0.0441238733 0.9779381
[2,] 0.0053055117 0.0106110234 0.9946945
[3,] 0.0010090105 0.0020180209 0.9989910
[4,] 0.0014124042 0.0028248083 0.9985876
[5,] 0.0004332930 0.0008665859 0.9995667
[6,] 0.0019903866 0.0039807732 0.9980096
[7,] 0.0513014203 0.1026028406 0.9486986
[8,] 0.1117893021 0.2235786042 0.8882107
[9,] 0.1215910501 0.2431821002 0.8784089
[10,] 0.1069410878 0.2138821756 0.8930589
[11,] 0.3574015462 0.7148030925 0.6425985
[12,] 0.2621141987 0.5242283974 0.7378858
[13,] 0.2764179456 0.5528358912 0.7235821
[14,] 0.1853528036 0.3707056073 0.8146472
[15,] 0.1377785820 0.2755571641 0.8622214
[16,] 0.2193708693 0.4387417387 0.7806291
[17,] 0.2751750952 0.5503501903 0.7248249
[18,] 0.3618615097 0.7237230194 0.6381385
[19,] 0.3190774042 0.6381548085 0.6809226
> postscript(file="/var/www/html/freestat/rcomp/tmp/1u66u1292963806.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/2u66u1292963806.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/3mfnx1292963806.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/4mfnx1292963806.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/5mfnx1292963806.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 = 60
Frequency = 1
1 2 3 4 5 6
-12.12669244 4.32766418 9.47605701 0.98167266 -8.49128109 2.05245790
7 8 9 10 11 12
-9.41468579 7.96206326 -2.79913166 -10.00328053 -1.84176045 0.24343776
13 14 15 16 17 18
2.80692136 6.58516546 3.25317804 4.84850414 -4.18248294 4.75837728
19 20 21 22 23 24
-5.04371633 6.43842258 4.66351245 -7.83992961 -0.12913402 -2.35791891
25 26 27 28 29 30
13.87178698 -1.66627625 -3.58110209 5.77817336 15.67733916 -5.19374585
31 32 33 34 35 36
13.13459041 -4.85382029 -0.09199703 0.47741475 1.70556381 7.14820750
37 38 39 40 41 42
-0.50595605 -10.37662424 -8.05156097 -6.61185278 -0.95751885 -13.09436526
43 44 45 46 47 48
-4.54963770 -5.51602290 -3.65495410 11.73929859 -4.17086625 -1.93099017
49 50 51 52 53 54
-4.04605985 1.13007085 -1.09657200 -4.99649738 -2.04605627 11.47727592
55 56 57 58 59 60
5.87344941 -4.03064265 1.88257035 5.62649680 4.43619690 -3.10273618
> postscript(file="/var/www/html/freestat/rcomp/tmp/6f7m01292963806.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -12.12669244 NA
1 4.32766418 -12.12669244
2 9.47605701 4.32766418
3 0.98167266 9.47605701
4 -8.49128109 0.98167266
5 2.05245790 -8.49128109
6 -9.41468579 2.05245790
7 7.96206326 -9.41468579
8 -2.79913166 7.96206326
9 -10.00328053 -2.79913166
10 -1.84176045 -10.00328053
11 0.24343776 -1.84176045
12 2.80692136 0.24343776
13 6.58516546 2.80692136
14 3.25317804 6.58516546
15 4.84850414 3.25317804
16 -4.18248294 4.84850414
17 4.75837728 -4.18248294
18 -5.04371633 4.75837728
19 6.43842258 -5.04371633
20 4.66351245 6.43842258
21 -7.83992961 4.66351245
22 -0.12913402 -7.83992961
23 -2.35791891 -0.12913402
24 13.87178698 -2.35791891
25 -1.66627625 13.87178698
26 -3.58110209 -1.66627625
27 5.77817336 -3.58110209
28 15.67733916 5.77817336
29 -5.19374585 15.67733916
30 13.13459041 -5.19374585
31 -4.85382029 13.13459041
32 -0.09199703 -4.85382029
33 0.47741475 -0.09199703
34 1.70556381 0.47741475
35 7.14820750 1.70556381
36 -0.50595605 7.14820750
37 -10.37662424 -0.50595605
38 -8.05156097 -10.37662424
39 -6.61185278 -8.05156097
40 -0.95751885 -6.61185278
41 -13.09436526 -0.95751885
42 -4.54963770 -13.09436526
43 -5.51602290 -4.54963770
44 -3.65495410 -5.51602290
45 11.73929859 -3.65495410
46 -4.17086625 11.73929859
47 -1.93099017 -4.17086625
48 -4.04605985 -1.93099017
49 1.13007085 -4.04605985
50 -1.09657200 1.13007085
51 -4.99649738 -1.09657200
52 -2.04605627 -4.99649738
53 11.47727592 -2.04605627
54 5.87344941 11.47727592
55 -4.03064265 5.87344941
56 1.88257035 -4.03064265
57 5.62649680 1.88257035
58 4.43619690 5.62649680
59 -3.10273618 4.43619690
60 NA -3.10273618
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 4.32766418 -12.12669244
[2,] 9.47605701 4.32766418
[3,] 0.98167266 9.47605701
[4,] -8.49128109 0.98167266
[5,] 2.05245790 -8.49128109
[6,] -9.41468579 2.05245790
[7,] 7.96206326 -9.41468579
[8,] -2.79913166 7.96206326
[9,] -10.00328053 -2.79913166
[10,] -1.84176045 -10.00328053
[11,] 0.24343776 -1.84176045
[12,] 2.80692136 0.24343776
[13,] 6.58516546 2.80692136
[14,] 3.25317804 6.58516546
[15,] 4.84850414 3.25317804
[16,] -4.18248294 4.84850414
[17,] 4.75837728 -4.18248294
[18,] -5.04371633 4.75837728
[19,] 6.43842258 -5.04371633
[20,] 4.66351245 6.43842258
[21,] -7.83992961 4.66351245
[22,] -0.12913402 -7.83992961
[23,] -2.35791891 -0.12913402
[24,] 13.87178698 -2.35791891
[25,] -1.66627625 13.87178698
[26,] -3.58110209 -1.66627625
[27,] 5.77817336 -3.58110209
[28,] 15.67733916 5.77817336
[29,] -5.19374585 15.67733916
[30,] 13.13459041 -5.19374585
[31,] -4.85382029 13.13459041
[32,] -0.09199703 -4.85382029
[33,] 0.47741475 -0.09199703
[34,] 1.70556381 0.47741475
[35,] 7.14820750 1.70556381
[36,] -0.50595605 7.14820750
[37,] -10.37662424 -0.50595605
[38,] -8.05156097 -10.37662424
[39,] -6.61185278 -8.05156097
[40,] -0.95751885 -6.61185278
[41,] -13.09436526 -0.95751885
[42,] -4.54963770 -13.09436526
[43,] -5.51602290 -4.54963770
[44,] -3.65495410 -5.51602290
[45,] 11.73929859 -3.65495410
[46,] -4.17086625 11.73929859
[47,] -1.93099017 -4.17086625
[48,] -4.04605985 -1.93099017
[49,] 1.13007085 -4.04605985
[50,] -1.09657200 1.13007085
[51,] -4.99649738 -1.09657200
[52,] -2.04605627 -4.99649738
[53,] 11.47727592 -2.04605627
[54,] 5.87344941 11.47727592
[55,] -4.03064265 5.87344941
[56,] 1.88257035 -4.03064265
[57,] 5.62649680 1.88257035
[58,] 4.43619690 5.62649680
[59,] -3.10273618 4.43619690
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 4.32766418 -12.12669244
2 9.47605701 4.32766418
3 0.98167266 9.47605701
4 -8.49128109 0.98167266
5 2.05245790 -8.49128109
6 -9.41468579 2.05245790
7 7.96206326 -9.41468579
8 -2.79913166 7.96206326
9 -10.00328053 -2.79913166
10 -1.84176045 -10.00328053
11 0.24343776 -1.84176045
12 2.80692136 0.24343776
13 6.58516546 2.80692136
14 3.25317804 6.58516546
15 4.84850414 3.25317804
16 -4.18248294 4.84850414
17 4.75837728 -4.18248294
18 -5.04371633 4.75837728
19 6.43842258 -5.04371633
20 4.66351245 6.43842258
21 -7.83992961 4.66351245
22 -0.12913402 -7.83992961
23 -2.35791891 -0.12913402
24 13.87178698 -2.35791891
25 -1.66627625 13.87178698
26 -3.58110209 -1.66627625
27 5.77817336 -3.58110209
28 15.67733916 5.77817336
29 -5.19374585 15.67733916
30 13.13459041 -5.19374585
31 -4.85382029 13.13459041
32 -0.09199703 -4.85382029
33 0.47741475 -0.09199703
34 1.70556381 0.47741475
35 7.14820750 1.70556381
36 -0.50595605 7.14820750
37 -10.37662424 -0.50595605
38 -8.05156097 -10.37662424
39 -6.61185278 -8.05156097
40 -0.95751885 -6.61185278
41 -13.09436526 -0.95751885
42 -4.54963770 -13.09436526
43 -5.51602290 -4.54963770
44 -3.65495410 -5.51602290
45 11.73929859 -3.65495410
46 -4.17086625 11.73929859
47 -1.93099017 -4.17086625
48 -4.04605985 -1.93099017
49 1.13007085 -4.04605985
50 -1.09657200 1.13007085
51 -4.99649738 -1.09657200
52 -2.04605627 -4.99649738
53 11.47727592 -2.04605627
54 5.87344941 11.47727592
55 -4.03064265 5.87344941
56 1.88257035 -4.03064265
57 5.62649680 1.88257035
58 4.43619690 5.62649680
59 -3.10273618 4.43619690
> 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/7qy3l1292963806.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/8qy3l1292963806.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/9qy3l1292963806.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/10i7361292963806.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/1148ju1292963806.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/127q0h1292963806.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/13lixq1292963806.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/1470ww1292963806.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/15a1u21292963806.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/16ekt81292963806.tab")
+ }
>
> try(system("convert tmp/1u66u1292963806.ps tmp/1u66u1292963806.png",intern=TRUE))
character(0)
> try(system("convert tmp/2u66u1292963806.ps tmp/2u66u1292963806.png",intern=TRUE))
character(0)
> try(system("convert tmp/3mfnx1292963806.ps tmp/3mfnx1292963806.png",intern=TRUE))
character(0)
> try(system("convert tmp/4mfnx1292963806.ps tmp/4mfnx1292963806.png",intern=TRUE))
character(0)
> try(system("convert tmp/5mfnx1292963806.ps tmp/5mfnx1292963806.png",intern=TRUE))
character(0)
> try(system("convert tmp/6f7m01292963806.ps tmp/6f7m01292963806.png",intern=TRUE))
character(0)
> try(system("convert tmp/7qy3l1292963806.ps tmp/7qy3l1292963806.png",intern=TRUE))
character(0)
> try(system("convert tmp/8qy3l1292963806.ps tmp/8qy3l1292963806.png",intern=TRUE))
character(0)
> try(system("convert tmp/9qy3l1292963806.ps tmp/9qy3l1292963806.png",intern=TRUE))
character(0)
> try(system("convert tmp/10i7361292963806.ps tmp/10i7361292963806.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.769 2.479 4.116