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(1.58,0.55,1.59,0.55,1.6,0.55,1.6,0.55,1.6,0.55,1.6,0.56,1.61,0.56,1.61,0.56,1.62,0.56,1.63,0.56,1.63,0.55,1.63,0.56,1.63,0.55,1.63,0.55,1.64,0.56,1.64,0.55,1.64,0.55,1.65,0.55,1.65,0.55,1.65,0.53,1.65,0.53,1.65,0.53,1.66,0.53,1.67,0.54,1.68,0.54,1.68,0.54,1.68,0.55,1.68,0.55,1.69,0.54,1.7,0.55,1.7,0.56,1.71,0.58,1.73,0.59,1.73,0.6,1.73,0.6,1.74,0.6,1.74,0.59,1.74,0.6,1.75,0.6,1.78,0.62,1.82,0.65,1.83,0.68,1.84,0.73,1.85,0.78,1.86,0.78,1.86,0.82,1.87,0.82,1.87,0.81,1.87,0.83,1.87,0.85,1.87,0.86,1.87,0.85,1.87,0.85,1.88,0.82,1.88,0.8,1.87,0.81,1.87,0.8,1.87,0.8,1.87,0.8,1.87,0.8,1.87,0.79),dim=c(2,61),dimnames=list(c('Y','X'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('Y','X'),1:61))
> 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.58 0.55 1 0 0 0 0 0 0 0 0 0 0 1
2 1.59 0.55 0 1 0 0 0 0 0 0 0 0 0 2
3 1.60 0.55 0 0 1 0 0 0 0 0 0 0 0 3
4 1.60 0.55 0 0 0 1 0 0 0 0 0 0 0 4
5 1.60 0.55 0 0 0 0 1 0 0 0 0 0 0 5
6 1.60 0.56 0 0 0 0 0 1 0 0 0 0 0 6
7 1.61 0.56 0 0 0 0 0 0 1 0 0 0 0 7
8 1.61 0.56 0 0 0 0 0 0 0 1 0 0 0 8
9 1.62 0.56 0 0 0 0 0 0 0 0 1 0 0 9
10 1.63 0.56 0 0 0 0 0 0 0 0 0 1 0 10
11 1.63 0.55 0 0 0 0 0 0 0 0 0 0 1 11
12 1.63 0.56 0 0 0 0 0 0 0 0 0 0 0 12
13 1.63 0.55 1 0 0 0 0 0 0 0 0 0 0 13
14 1.63 0.55 0 1 0 0 0 0 0 0 0 0 0 14
15 1.64 0.56 0 0 1 0 0 0 0 0 0 0 0 15
16 1.64 0.55 0 0 0 1 0 0 0 0 0 0 0 16
17 1.64 0.55 0 0 0 0 1 0 0 0 0 0 0 17
18 1.65 0.55 0 0 0 0 0 1 0 0 0 0 0 18
19 1.65 0.55 0 0 0 0 0 0 1 0 0 0 0 19
20 1.65 0.53 0 0 0 0 0 0 0 1 0 0 0 20
21 1.65 0.53 0 0 0 0 0 0 0 0 1 0 0 21
22 1.65 0.53 0 0 0 0 0 0 0 0 0 1 0 22
23 1.66 0.53 0 0 0 0 0 0 0 0 0 0 1 23
24 1.67 0.54 0 0 0 0 0 0 0 0 0 0 0 24
25 1.68 0.54 1 0 0 0 0 0 0 0 0 0 0 25
26 1.68 0.54 0 1 0 0 0 0 0 0 0 0 0 26
27 1.68 0.55 0 0 1 0 0 0 0 0 0 0 0 27
28 1.68 0.55 0 0 0 1 0 0 0 0 0 0 0 28
29 1.69 0.54 0 0 0 0 1 0 0 0 0 0 0 29
30 1.70 0.55 0 0 0 0 0 1 0 0 0 0 0 30
31 1.70 0.56 0 0 0 0 0 0 1 0 0 0 0 31
32 1.71 0.58 0 0 0 0 0 0 0 1 0 0 0 32
33 1.73 0.59 0 0 0 0 0 0 0 0 1 0 0 33
34 1.73 0.60 0 0 0 0 0 0 0 0 0 1 0 34
35 1.73 0.60 0 0 0 0 0 0 0 0 0 0 1 35
36 1.74 0.60 0 0 0 0 0 0 0 0 0 0 0 36
37 1.74 0.59 1 0 0 0 0 0 0 0 0 0 0 37
38 1.74 0.60 0 1 0 0 0 0 0 0 0 0 0 38
39 1.75 0.60 0 0 1 0 0 0 0 0 0 0 0 39
40 1.78 0.62 0 0 0 1 0 0 0 0 0 0 0 40
41 1.82 0.65 0 0 0 0 1 0 0 0 0 0 0 41
42 1.83 0.68 0 0 0 0 0 1 0 0 0 0 0 42
43 1.84 0.73 0 0 0 0 0 0 1 0 0 0 0 43
44 1.85 0.78 0 0 0 0 0 0 0 1 0 0 0 44
45 1.86 0.78 0 0 0 0 0 0 0 0 1 0 0 45
46 1.86 0.82 0 0 0 0 0 0 0 0 0 1 0 46
47 1.87 0.82 0 0 0 0 0 0 0 0 0 0 1 47
48 1.87 0.81 0 0 0 0 0 0 0 0 0 0 0 48
49 1.87 0.83 1 0 0 0 0 0 0 0 0 0 0 49
50 1.87 0.85 0 1 0 0 0 0 0 0 0 0 0 50
51 1.87 0.86 0 0 1 0 0 0 0 0 0 0 0 51
52 1.87 0.85 0 0 0 1 0 0 0 0 0 0 0 52
53 1.87 0.85 0 0 0 0 1 0 0 0 0 0 0 53
54 1.88 0.82 0 0 0 0 0 1 0 0 0 0 0 54
55 1.88 0.80 0 0 0 0 0 0 1 0 0 0 0 55
56 1.87 0.81 0 0 0 0 0 0 0 1 0 0 0 56
57 1.87 0.80 0 0 0 0 0 0 0 0 1 0 0 57
58 1.87 0.80 0 0 0 0 0 0 0 0 0 1 0 58
59 1.87 0.80 0 0 0 0 0 0 0 0 0 0 1 59
60 1.87 0.80 0 0 0 0 0 0 0 0 0 0 0 60
61 1.87 0.79 1 0 0 0 0 0 0 0 0 0 0 61
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
1.3945018 0.3401696 -0.0018185 -0.0011698 -0.0009971 0.0012167
M5 M6 M7 M8 M9 M10
0.0060697 0.0089227 0.0064151 0.0005468 0.0047605 -0.0004274
M11 t
0.0004666 0.0037863
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.023814 -0.008671 -0.001790 0.005091 0.043081
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.3945018 0.0166442 83.783 < 2e-16 ***
X 0.3401696 0.0322041 10.563 5.31e-14 ***
M1 -0.0018185 0.0092558 -0.196 0.845
M2 -0.0011698 0.0097206 -0.120 0.905
M3 -0.0009971 0.0097086 -0.103 0.919
M4 0.0012167 0.0096907 0.126 0.901
M5 0.0060697 0.0096795 0.627 0.534
M6 0.0089227 0.0096700 0.923 0.361
M7 0.0064151 0.0096641 0.664 0.510
M8 0.0005468 0.0096649 0.057 0.955
M9 0.0047605 0.0096540 0.493 0.624
M10 -0.0004274 0.0096550 -0.044 0.965
M11 0.0004666 0.0096466 0.048 0.962
t 0.0037863 0.0002172 17.433 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.01525 on 47 degrees of freedom
Multiple R-squared: 0.9834, Adjusted R-squared: 0.9789
F-statistic: 214.8 on 13 and 47 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,] 1.075725e-02 2.151449e-02 0.98924275
[2,] 2.328762e-03 4.657524e-03 0.99767124
[3,] 1.115960e-03 2.231919e-03 0.99888404
[4,] 3.447160e-04 6.894320e-04 0.99965528
[5,] 3.519455e-04 7.038910e-04 0.99964805
[6,] 8.571142e-04 1.714228e-03 0.99914289
[7,] 3.199415e-04 6.398831e-04 0.99968006
[8,] 1.180628e-04 2.361256e-04 0.99988194
[9,] 2.004077e-04 4.008154e-04 0.99979959
[10,] 8.087453e-05 1.617491e-04 0.99991913
[11,] 3.474447e-05 6.948893e-05 0.99996526
[12,] 2.180481e-05 4.360961e-05 0.99997820
[13,] 1.543151e-05 3.086302e-05 0.99998457
[14,] 2.272023e-05 4.544047e-05 0.99997728
[15,] 3.122569e-05 6.245139e-05 0.99996877
[16,] 2.675149e-05 5.350297e-05 0.99997325
[17,] 3.229067e-05 6.458133e-05 0.99996771
[18,] 2.424719e-05 4.849438e-05 0.99997575
[19,] 8.571389e-05 1.714278e-04 0.99991429
[20,] 2.326964e-04 4.653928e-04 0.99976730
[21,] 1.401922e-03 2.803843e-03 0.99859808
[22,] 1.314155e-02 2.628309e-02 0.98685845
[23,] 1.303732e-01 2.607464e-01 0.86962682
[24,] 5.755731e-01 8.488539e-01 0.42442694
[25,] 9.387233e-01 1.225534e-01 0.06127672
[26,] 9.190959e-01 1.618082e-01 0.08090411
[27,] 9.245293e-01 1.509413e-01 0.07547066
[28,] 9.162073e-01 1.675854e-01 0.08379270
> postscript(file="/var/www/html/rcomp/tmp/1x3n91258717790.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/rcomp/tmp/25heg1258717790.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/rcomp/tmp/31v961258717790.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/rcomp/tmp/4l03a1258717790.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/rcomp/tmp/5w09d1258717790.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 = 61
Frequency = 1
1 2 3 4 5
-0.0035628224 0.0020021707 0.0080431885 0.0020431885 -0.0065961330
6 7 8 9 10
-0.0166371508 -0.0079157938 -0.0058337583 -0.0038337583 0.0075679380
11 12 13 14 15
0.0062892950 -0.0004320620 0.0010018595 -0.0034331474 -0.0007938259
16 17 18 19 20
-0.0033921296 -0.0120314511 -0.0086707726 -0.0099494156 -0.0010639876
21 22 23 24 25
-0.0090639876 -0.0076622914 -0.0023426306 0.0009360124 0.0089682376
26 27 28 29 30
0.0045332307 -0.0028274478 -0.0088274478 -0.0040650730 -0.0041060908
31 32 33 34 35
-0.0087864300 -0.0035077870 0.0050905167 0.0030905167 -0.0015898225
36 37 38 39 40
0.0050905167 0.0065244382 -0.0013122649 0.0047287529 0.0219253604
41 42 43 44 45
0.0430809501 0.0362365399 0.0279494156 0.0230229699 0.0250229699
46 47 48 49 50
0.0128178811 0.0181375419 0.0182195774 0.0094484101 -0.0017899892
51 52 53 54 55
-0.0091506677 -0.0117489715 -0.0203882930 -0.0068225257 -0.0012977762
56 57 58 59 60
-0.0126174370 -0.0172157407 -0.0158140445 -0.0204943837 -0.0238140445
61
-0.0223801230
> postscript(file="/var/www/html/rcomp/tmp/6sijo1258717790.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.0035628224 NA
1 0.0020021707 -0.0035628224
2 0.0080431885 0.0020021707
3 0.0020431885 0.0080431885
4 -0.0065961330 0.0020431885
5 -0.0166371508 -0.0065961330
6 -0.0079157938 -0.0166371508
7 -0.0058337583 -0.0079157938
8 -0.0038337583 -0.0058337583
9 0.0075679380 -0.0038337583
10 0.0062892950 0.0075679380
11 -0.0004320620 0.0062892950
12 0.0010018595 -0.0004320620
13 -0.0034331474 0.0010018595
14 -0.0007938259 -0.0034331474
15 -0.0033921296 -0.0007938259
16 -0.0120314511 -0.0033921296
17 -0.0086707726 -0.0120314511
18 -0.0099494156 -0.0086707726
19 -0.0010639876 -0.0099494156
20 -0.0090639876 -0.0010639876
21 -0.0076622914 -0.0090639876
22 -0.0023426306 -0.0076622914
23 0.0009360124 -0.0023426306
24 0.0089682376 0.0009360124
25 0.0045332307 0.0089682376
26 -0.0028274478 0.0045332307
27 -0.0088274478 -0.0028274478
28 -0.0040650730 -0.0088274478
29 -0.0041060908 -0.0040650730
30 -0.0087864300 -0.0041060908
31 -0.0035077870 -0.0087864300
32 0.0050905167 -0.0035077870
33 0.0030905167 0.0050905167
34 -0.0015898225 0.0030905167
35 0.0050905167 -0.0015898225
36 0.0065244382 0.0050905167
37 -0.0013122649 0.0065244382
38 0.0047287529 -0.0013122649
39 0.0219253604 0.0047287529
40 0.0430809501 0.0219253604
41 0.0362365399 0.0430809501
42 0.0279494156 0.0362365399
43 0.0230229699 0.0279494156
44 0.0250229699 0.0230229699
45 0.0128178811 0.0250229699
46 0.0181375419 0.0128178811
47 0.0182195774 0.0181375419
48 0.0094484101 0.0182195774
49 -0.0017899892 0.0094484101
50 -0.0091506677 -0.0017899892
51 -0.0117489715 -0.0091506677
52 -0.0203882930 -0.0117489715
53 -0.0068225257 -0.0203882930
54 -0.0012977762 -0.0068225257
55 -0.0126174370 -0.0012977762
56 -0.0172157407 -0.0126174370
57 -0.0158140445 -0.0172157407
58 -0.0204943837 -0.0158140445
59 -0.0238140445 -0.0204943837
60 -0.0223801230 -0.0238140445
61 NA -0.0223801230
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.0020021707 -0.0035628224
[2,] 0.0080431885 0.0020021707
[3,] 0.0020431885 0.0080431885
[4,] -0.0065961330 0.0020431885
[5,] -0.0166371508 -0.0065961330
[6,] -0.0079157938 -0.0166371508
[7,] -0.0058337583 -0.0079157938
[8,] -0.0038337583 -0.0058337583
[9,] 0.0075679380 -0.0038337583
[10,] 0.0062892950 0.0075679380
[11,] -0.0004320620 0.0062892950
[12,] 0.0010018595 -0.0004320620
[13,] -0.0034331474 0.0010018595
[14,] -0.0007938259 -0.0034331474
[15,] -0.0033921296 -0.0007938259
[16,] -0.0120314511 -0.0033921296
[17,] -0.0086707726 -0.0120314511
[18,] -0.0099494156 -0.0086707726
[19,] -0.0010639876 -0.0099494156
[20,] -0.0090639876 -0.0010639876
[21,] -0.0076622914 -0.0090639876
[22,] -0.0023426306 -0.0076622914
[23,] 0.0009360124 -0.0023426306
[24,] 0.0089682376 0.0009360124
[25,] 0.0045332307 0.0089682376
[26,] -0.0028274478 0.0045332307
[27,] -0.0088274478 -0.0028274478
[28,] -0.0040650730 -0.0088274478
[29,] -0.0041060908 -0.0040650730
[30,] -0.0087864300 -0.0041060908
[31,] -0.0035077870 -0.0087864300
[32,] 0.0050905167 -0.0035077870
[33,] 0.0030905167 0.0050905167
[34,] -0.0015898225 0.0030905167
[35,] 0.0050905167 -0.0015898225
[36,] 0.0065244382 0.0050905167
[37,] -0.0013122649 0.0065244382
[38,] 0.0047287529 -0.0013122649
[39,] 0.0219253604 0.0047287529
[40,] 0.0430809501 0.0219253604
[41,] 0.0362365399 0.0430809501
[42,] 0.0279494156 0.0362365399
[43,] 0.0230229699 0.0279494156
[44,] 0.0250229699 0.0230229699
[45,] 0.0128178811 0.0250229699
[46,] 0.0181375419 0.0128178811
[47,] 0.0182195774 0.0181375419
[48,] 0.0094484101 0.0182195774
[49,] -0.0017899892 0.0094484101
[50,] -0.0091506677 -0.0017899892
[51,] -0.0117489715 -0.0091506677
[52,] -0.0203882930 -0.0117489715
[53,] -0.0068225257 -0.0203882930
[54,] -0.0012977762 -0.0068225257
[55,] -0.0126174370 -0.0012977762
[56,] -0.0172157407 -0.0126174370
[57,] -0.0158140445 -0.0172157407
[58,] -0.0204943837 -0.0158140445
[59,] -0.0238140445 -0.0204943837
[60,] -0.0223801230 -0.0238140445
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.0020021707 -0.0035628224
2 0.0080431885 0.0020021707
3 0.0020431885 0.0080431885
4 -0.0065961330 0.0020431885
5 -0.0166371508 -0.0065961330
6 -0.0079157938 -0.0166371508
7 -0.0058337583 -0.0079157938
8 -0.0038337583 -0.0058337583
9 0.0075679380 -0.0038337583
10 0.0062892950 0.0075679380
11 -0.0004320620 0.0062892950
12 0.0010018595 -0.0004320620
13 -0.0034331474 0.0010018595
14 -0.0007938259 -0.0034331474
15 -0.0033921296 -0.0007938259
16 -0.0120314511 -0.0033921296
17 -0.0086707726 -0.0120314511
18 -0.0099494156 -0.0086707726
19 -0.0010639876 -0.0099494156
20 -0.0090639876 -0.0010639876
21 -0.0076622914 -0.0090639876
22 -0.0023426306 -0.0076622914
23 0.0009360124 -0.0023426306
24 0.0089682376 0.0009360124
25 0.0045332307 0.0089682376
26 -0.0028274478 0.0045332307
27 -0.0088274478 -0.0028274478
28 -0.0040650730 -0.0088274478
29 -0.0041060908 -0.0040650730
30 -0.0087864300 -0.0041060908
31 -0.0035077870 -0.0087864300
32 0.0050905167 -0.0035077870
33 0.0030905167 0.0050905167
34 -0.0015898225 0.0030905167
35 0.0050905167 -0.0015898225
36 0.0065244382 0.0050905167
37 -0.0013122649 0.0065244382
38 0.0047287529 -0.0013122649
39 0.0219253604 0.0047287529
40 0.0430809501 0.0219253604
41 0.0362365399 0.0430809501
42 0.0279494156 0.0362365399
43 0.0230229699 0.0279494156
44 0.0250229699 0.0230229699
45 0.0128178811 0.0250229699
46 0.0181375419 0.0128178811
47 0.0182195774 0.0181375419
48 0.0094484101 0.0182195774
49 -0.0017899892 0.0094484101
50 -0.0091506677 -0.0017899892
51 -0.0117489715 -0.0091506677
52 -0.0203882930 -0.0117489715
53 -0.0068225257 -0.0203882930
54 -0.0012977762 -0.0068225257
55 -0.0126174370 -0.0012977762
56 -0.0172157407 -0.0126174370
57 -0.0158140445 -0.0172157407
58 -0.0204943837 -0.0158140445
59 -0.0238140445 -0.0204943837
60 -0.0223801230 -0.0238140445
> 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/78urh1258717791.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/rcomp/tmp/8vlji1258717791.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/rcomp/tmp/90dtk1258717791.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/rcomp/tmp/101nol1258717791.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/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/111vor1258717791.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/12clpz1258717791.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/13rdve1258717791.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/14y7cz1258717791.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/15dn331258717791.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/166anp1258717791.tab")
+ }
>
> system("convert tmp/1x3n91258717790.ps tmp/1x3n91258717790.png")
> system("convert tmp/25heg1258717790.ps tmp/25heg1258717790.png")
> system("convert tmp/31v961258717790.ps tmp/31v961258717790.png")
> system("convert tmp/4l03a1258717790.ps tmp/4l03a1258717790.png")
> system("convert tmp/5w09d1258717790.ps tmp/5w09d1258717790.png")
> system("convert tmp/6sijo1258717790.ps tmp/6sijo1258717790.png")
> system("convert tmp/78urh1258717791.ps tmp/78urh1258717791.png")
> system("convert tmp/8vlji1258717791.ps tmp/8vlji1258717791.png")
> system("convert tmp/90dtk1258717791.ps tmp/90dtk1258717791.png")
> system("convert tmp/101nol1258717791.ps tmp/101nol1258717791.png")
>
>
> proc.time()
user system elapsed
2.405 1.563 2.799