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(7.2,-6,7.5,8.3,7.4,0,7.2,7.5,8.8,-4,7.4,7.2,9.3,-2,8.8,7.4,9.3,-2,9.3,8.8,8.7,-6,9.3,9.3,8.2,-7,8.7,9.3,8.3,-6,8.2,8.7,8.5,-6,8.3,8.2,8.6,-3,8.5,8.3,8.5,-2,8.6,8.5,8.2,-5,8.5,8.6,8.1,-11,8.2,8.5,7.9,-11,8.1,8.2,8.6,-11,7.9,8.1,8.7,-10,8.6,7.9,8.7,-14,8.7,8.6,8.5,-8,8.7,8.7,8.4,-9,8.5,8.7,8.5,-5,8.4,8.5,8.7,-1,8.5,8.4,8.7,-2,8.7,8.5,8.6,-5,8.7,8.7,8.5,-4,8.6,8.7,8.3,-6,8.5,8.6,8,-2,8.3,8.5,8.2,-2,8,8.3,8.1,-2,8.2,8,8.1,-2,8.1,8.2,8,2,8.1,8.1,7.9,1,8,8.1,7.9,-8,7.9,8,8,-1,7.9,7.9,8,1,8,7.9,7.9,-1,8,8,8,2,7.9,8,7.7,2,8,7.9,7.2,1,7.7,8,7.5,-1,7.2,7.7,7.3,-2,7.5,7.2,7,-2,7.3,7.5,7,-1,7,7.3,7,-8,7,7,7.2,-4,7,7,7.3,-6,7.2,7,7.1,-3,7.3,7.2,6.8,-3,7.1,7.3,6.4,-7,6.8,7.1,6.1,-9,6.4,6.8,6.5,-11,6.1,6.4,7.7,-13,6.5,6.1,7.9,-11,7.7,6.5,7.5,-9,7.9,7.7,6.9,-17,7.5,7.9,6.6,-22,6.9,7.5,6.9,-25,6.6,6.9),dim=c(4,56),dimnames=list(c('TW','CV','TW1','TW2'),1:56))
> y <- array(NA,dim=c(4,56),dimnames=list(c('TW','CV','TW1','TW2'),1:56))
> 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
TW CV TW1 TW2 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 7.2 -6 7.5 8.3 1 0 0 0 0 0 0 0 0 0 0 1
2 7.4 0 7.2 7.5 0 1 0 0 0 0 0 0 0 0 0 2
3 8.8 -4 7.4 7.2 0 0 1 0 0 0 0 0 0 0 0 3
4 9.3 -2 8.8 7.4 0 0 0 1 0 0 0 0 0 0 0 4
5 9.3 -2 9.3 8.8 0 0 0 0 1 0 0 0 0 0 0 5
6 8.7 -6 9.3 9.3 0 0 0 0 0 1 0 0 0 0 0 6
7 8.2 -7 8.7 9.3 0 0 0 0 0 0 1 0 0 0 0 7
8 8.3 -6 8.2 8.7 0 0 0 0 0 0 0 1 0 0 0 8
9 8.5 -6 8.3 8.2 0 0 0 0 0 0 0 0 1 0 0 9
10 8.6 -3 8.5 8.3 0 0 0 0 0 0 0 0 0 1 0 10
11 8.5 -2 8.6 8.5 0 0 0 0 0 0 0 0 0 0 1 11
12 8.2 -5 8.5 8.6 0 0 0 0 0 0 0 0 0 0 0 12
13 8.1 -11 8.2 8.5 1 0 0 0 0 0 0 0 0 0 0 13
14 7.9 -11 8.1 8.2 0 1 0 0 0 0 0 0 0 0 0 14
15 8.6 -11 7.9 8.1 0 0 1 0 0 0 0 0 0 0 0 15
16 8.7 -10 8.6 7.9 0 0 0 1 0 0 0 0 0 0 0 16
17 8.7 -14 8.7 8.6 0 0 0 0 1 0 0 0 0 0 0 17
18 8.5 -8 8.7 8.7 0 0 0 0 0 1 0 0 0 0 0 18
19 8.4 -9 8.5 8.7 0 0 0 0 0 0 1 0 0 0 0 19
20 8.5 -5 8.4 8.5 0 0 0 0 0 0 0 1 0 0 0 20
21 8.7 -1 8.5 8.4 0 0 0 0 0 0 0 0 1 0 0 21
22 8.7 -2 8.7 8.5 0 0 0 0 0 0 0 0 0 1 0 22
23 8.6 -5 8.7 8.7 0 0 0 0 0 0 0 0 0 0 1 23
24 8.5 -4 8.6 8.7 0 0 0 0 0 0 0 0 0 0 0 24
25 8.3 -6 8.5 8.6 1 0 0 0 0 0 0 0 0 0 0 25
26 8.0 -2 8.3 8.5 0 1 0 0 0 0 0 0 0 0 0 26
27 8.2 -2 8.0 8.3 0 0 1 0 0 0 0 0 0 0 0 27
28 8.1 -2 8.2 8.0 0 0 0 1 0 0 0 0 0 0 0 28
29 8.1 -2 8.1 8.2 0 0 0 0 1 0 0 0 0 0 0 29
30 8.0 2 8.1 8.1 0 0 0 0 0 1 0 0 0 0 0 30
31 7.9 1 8.0 8.1 0 0 0 0 0 0 1 0 0 0 0 31
32 7.9 -8 7.9 8.0 0 0 0 0 0 0 0 1 0 0 0 32
33 8.0 -1 7.9 7.9 0 0 0 0 0 0 0 0 1 0 0 33
34 8.0 1 8.0 7.9 0 0 0 0 0 0 0 0 0 1 0 34
35 7.9 -1 8.0 8.0 0 0 0 0 0 0 0 0 0 0 1 35
36 8.0 2 7.9 8.0 0 0 0 0 0 0 0 0 0 0 0 36
37 7.7 2 8.0 7.9 1 0 0 0 0 0 0 0 0 0 0 37
38 7.2 1 7.7 8.0 0 1 0 0 0 0 0 0 0 0 0 38
39 7.5 -1 7.2 7.7 0 0 1 0 0 0 0 0 0 0 0 39
40 7.3 -2 7.5 7.2 0 0 0 1 0 0 0 0 0 0 0 40
41 7.0 -2 7.3 7.5 0 0 0 0 1 0 0 0 0 0 0 41
42 7.0 -1 7.0 7.3 0 0 0 0 0 1 0 0 0 0 0 42
43 7.0 -8 7.0 7.0 0 0 0 0 0 0 1 0 0 0 0 43
44 7.2 -4 7.0 7.0 0 0 0 0 0 0 0 1 0 0 0 44
45 7.3 -6 7.2 7.0 0 0 0 0 0 0 0 0 1 0 0 45
46 7.1 -3 7.3 7.2 0 0 0 0 0 0 0 0 0 1 0 46
47 6.8 -3 7.1 7.3 0 0 0 0 0 0 0 0 0 0 1 47
48 6.4 -7 6.8 7.1 0 0 0 0 0 0 0 0 0 0 0 48
49 6.1 -9 6.4 6.8 1 0 0 0 0 0 0 0 0 0 0 49
50 6.5 -11 6.1 6.4 0 1 0 0 0 0 0 0 0 0 0 50
51 7.7 -13 6.5 6.1 0 0 1 0 0 0 0 0 0 0 0 51
52 7.9 -11 7.7 6.5 0 0 0 1 0 0 0 0 0 0 0 52
53 7.5 -9 7.9 7.7 0 0 0 0 1 0 0 0 0 0 0 53
54 6.9 -17 7.5 7.9 0 0 0 0 0 1 0 0 0 0 0 54
55 6.6 -22 6.9 7.5 0 0 0 0 0 0 1 0 0 0 0 55
56 6.9 -25 6.6 6.9 0 0 0 0 0 0 0 1 0 0 0 56
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) CV TW1 TW2 M1 M2
2.757148 -0.003377 1.326237 -0.639603 -0.109036 -0.046001
M3 M4 M5 M6 M7 M8
0.673083 -0.271432 -0.047416 -0.086565 0.023523 0.246756
M9 M10 M11 t
0.137257 -0.004917 -0.017305 -0.011892
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-0.305181 -0.104774 0.002774 0.108968 0.349729
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.757148 0.527471 5.227 5.73e-06 ***
CV -0.003377 0.004699 -0.719 0.4766
TW1 1.326237 0.101230 13.101 4.68e-16 ***
TW2 -0.639603 0.099888 -6.403 1.28e-07 ***
M1 -0.109036 0.118356 -0.921 0.3624
M2 -0.046001 0.120025 -0.383 0.7036
M3 0.673083 0.122248 5.506 2.34e-06 ***
M4 -0.271432 0.145796 -1.862 0.0700 .
M5 -0.047416 0.118628 -0.400 0.6915
M6 -0.086565 0.116184 -0.745 0.4606
M7 0.023523 0.118638 0.198 0.8438
M8 0.246756 0.119056 2.073 0.0447 *
M9 0.137257 0.124156 1.106 0.2755
M10 -0.004917 0.124872 -0.039 0.9688
M11 -0.017305 0.121992 -0.142 0.8879
t -0.011892 0.002358 -5.043 1.03e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1716 on 40 degrees of freedom
Multiple R-squared: 0.9613, Adjusted R-squared: 0.9468
F-statistic: 66.24 on 15 and 40 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,] 0.17083671 0.34167341 0.8291633
[2,] 0.20151669 0.40303338 0.7984833
[3,] 0.15005090 0.30010180 0.8499491
[4,] 0.07592135 0.15184270 0.9240786
[5,] 0.04924111 0.09848222 0.9507589
[6,] 0.05157328 0.10314655 0.9484267
[7,] 0.02721143 0.05442286 0.9727886
[8,] 0.01492082 0.02984165 0.9850792
[9,] 0.15361242 0.30722483 0.8463876
[10,] 0.10273999 0.20547998 0.8972600
[11,] 0.07267349 0.14534699 0.9273265
[12,] 0.04676581 0.09353162 0.9532342
[13,] 0.04268574 0.08537149 0.9573143
[14,] 0.15218622 0.30437244 0.8478138
[15,] 0.10002098 0.20004197 0.8999790
[16,] 0.05904340 0.11808680 0.9409566
[17,] 0.03411727 0.06823455 0.9658827
[18,] 0.24041438 0.48082876 0.7595856
[19,] 0.75241724 0.49516553 0.2475828
> postscript(file="/var/www/html/rcomp/tmp/1xo2f1260821748.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/20y4x1260821748.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/3jdjl1260821748.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/4s5qi1260821748.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/5doy51260821748.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 = 56
Frequency = 1
1 2 3 4 5 6
-0.094548331 -0.039243121 0.182930384 -0.082720240 -0.062518176 -0.305181329
7 8 9 10 11 12
-0.111012499 0.060380210 -0.070654558 -0.007745567 -0.084791672 -0.203750593
13 14 15 16 17 18
0.130827928 -0.179572553 0.014522988 0.018020132 0.107487942 0.042750110
19 20 21 22 23 24
0.106424223 0.013293530 0.151606743 0.101009076 0.143079993 0.173667398
25 26 27 28 29 30
0.156505219 0.020155736 -0.217085374 0.182193171 0.230613344 0.131201522
31 32 33 34 35 36
0.062251955 -0.110815044 0.070251841 0.098447516 0.079934767 0.317275501
37 38 39 40 41 42
-0.058380708 -0.151069172 -0.093776583 -0.058419043 -0.113414860 0.210954032
43 44 45 46 47 48
-0.102760193 -0.100593903 -0.151204026 -0.191711025 -0.138223087 -0.287192306
49 50 51 52 53 54
-0.134404109 0.349729109 0.113408584 -0.059074020 -0.162168250 -0.079724336
55 56
0.045096514 0.137735206
> postscript(file="/var/www/html/rcomp/tmp/62of21260821748.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 = 56
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.094548331 NA
1 -0.039243121 -0.094548331
2 0.182930384 -0.039243121
3 -0.082720240 0.182930384
4 -0.062518176 -0.082720240
5 -0.305181329 -0.062518176
6 -0.111012499 -0.305181329
7 0.060380210 -0.111012499
8 -0.070654558 0.060380210
9 -0.007745567 -0.070654558
10 -0.084791672 -0.007745567
11 -0.203750593 -0.084791672
12 0.130827928 -0.203750593
13 -0.179572553 0.130827928
14 0.014522988 -0.179572553
15 0.018020132 0.014522988
16 0.107487942 0.018020132
17 0.042750110 0.107487942
18 0.106424223 0.042750110
19 0.013293530 0.106424223
20 0.151606743 0.013293530
21 0.101009076 0.151606743
22 0.143079993 0.101009076
23 0.173667398 0.143079993
24 0.156505219 0.173667398
25 0.020155736 0.156505219
26 -0.217085374 0.020155736
27 0.182193171 -0.217085374
28 0.230613344 0.182193171
29 0.131201522 0.230613344
30 0.062251955 0.131201522
31 -0.110815044 0.062251955
32 0.070251841 -0.110815044
33 0.098447516 0.070251841
34 0.079934767 0.098447516
35 0.317275501 0.079934767
36 -0.058380708 0.317275501
37 -0.151069172 -0.058380708
38 -0.093776583 -0.151069172
39 -0.058419043 -0.093776583
40 -0.113414860 -0.058419043
41 0.210954032 -0.113414860
42 -0.102760193 0.210954032
43 -0.100593903 -0.102760193
44 -0.151204026 -0.100593903
45 -0.191711025 -0.151204026
46 -0.138223087 -0.191711025
47 -0.287192306 -0.138223087
48 -0.134404109 -0.287192306
49 0.349729109 -0.134404109
50 0.113408584 0.349729109
51 -0.059074020 0.113408584
52 -0.162168250 -0.059074020
53 -0.079724336 -0.162168250
54 0.045096514 -0.079724336
55 0.137735206 0.045096514
56 NA 0.137735206
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.039243121 -0.094548331
[2,] 0.182930384 -0.039243121
[3,] -0.082720240 0.182930384
[4,] -0.062518176 -0.082720240
[5,] -0.305181329 -0.062518176
[6,] -0.111012499 -0.305181329
[7,] 0.060380210 -0.111012499
[8,] -0.070654558 0.060380210
[9,] -0.007745567 -0.070654558
[10,] -0.084791672 -0.007745567
[11,] -0.203750593 -0.084791672
[12,] 0.130827928 -0.203750593
[13,] -0.179572553 0.130827928
[14,] 0.014522988 -0.179572553
[15,] 0.018020132 0.014522988
[16,] 0.107487942 0.018020132
[17,] 0.042750110 0.107487942
[18,] 0.106424223 0.042750110
[19,] 0.013293530 0.106424223
[20,] 0.151606743 0.013293530
[21,] 0.101009076 0.151606743
[22,] 0.143079993 0.101009076
[23,] 0.173667398 0.143079993
[24,] 0.156505219 0.173667398
[25,] 0.020155736 0.156505219
[26,] -0.217085374 0.020155736
[27,] 0.182193171 -0.217085374
[28,] 0.230613344 0.182193171
[29,] 0.131201522 0.230613344
[30,] 0.062251955 0.131201522
[31,] -0.110815044 0.062251955
[32,] 0.070251841 -0.110815044
[33,] 0.098447516 0.070251841
[34,] 0.079934767 0.098447516
[35,] 0.317275501 0.079934767
[36,] -0.058380708 0.317275501
[37,] -0.151069172 -0.058380708
[38,] -0.093776583 -0.151069172
[39,] -0.058419043 -0.093776583
[40,] -0.113414860 -0.058419043
[41,] 0.210954032 -0.113414860
[42,] -0.102760193 0.210954032
[43,] -0.100593903 -0.102760193
[44,] -0.151204026 -0.100593903
[45,] -0.191711025 -0.151204026
[46,] -0.138223087 -0.191711025
[47,] -0.287192306 -0.138223087
[48,] -0.134404109 -0.287192306
[49,] 0.349729109 -0.134404109
[50,] 0.113408584 0.349729109
[51,] -0.059074020 0.113408584
[52,] -0.162168250 -0.059074020
[53,] -0.079724336 -0.162168250
[54,] 0.045096514 -0.079724336
[55,] 0.137735206 0.045096514
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.039243121 -0.094548331
2 0.182930384 -0.039243121
3 -0.082720240 0.182930384
4 -0.062518176 -0.082720240
5 -0.305181329 -0.062518176
6 -0.111012499 -0.305181329
7 0.060380210 -0.111012499
8 -0.070654558 0.060380210
9 -0.007745567 -0.070654558
10 -0.084791672 -0.007745567
11 -0.203750593 -0.084791672
12 0.130827928 -0.203750593
13 -0.179572553 0.130827928
14 0.014522988 -0.179572553
15 0.018020132 0.014522988
16 0.107487942 0.018020132
17 0.042750110 0.107487942
18 0.106424223 0.042750110
19 0.013293530 0.106424223
20 0.151606743 0.013293530
21 0.101009076 0.151606743
22 0.143079993 0.101009076
23 0.173667398 0.143079993
24 0.156505219 0.173667398
25 0.020155736 0.156505219
26 -0.217085374 0.020155736
27 0.182193171 -0.217085374
28 0.230613344 0.182193171
29 0.131201522 0.230613344
30 0.062251955 0.131201522
31 -0.110815044 0.062251955
32 0.070251841 -0.110815044
33 0.098447516 0.070251841
34 0.079934767 0.098447516
35 0.317275501 0.079934767
36 -0.058380708 0.317275501
37 -0.151069172 -0.058380708
38 -0.093776583 -0.151069172
39 -0.058419043 -0.093776583
40 -0.113414860 -0.058419043
41 0.210954032 -0.113414860
42 -0.102760193 0.210954032
43 -0.100593903 -0.102760193
44 -0.151204026 -0.100593903
45 -0.191711025 -0.151204026
46 -0.138223087 -0.191711025
47 -0.287192306 -0.138223087
48 -0.134404109 -0.287192306
49 0.349729109 -0.134404109
50 0.113408584 0.349729109
51 -0.059074020 0.113408584
52 -0.162168250 -0.059074020
53 -0.079724336 -0.162168250
54 0.045096514 -0.079724336
55 0.137735206 0.045096514
> 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/7la1e1260821748.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/8a0ml1260821748.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/9wpgc1260821748.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/106hd31260821748.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/11f29l1260821748.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/12w2sg1260821748.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/138uee1260821748.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/140i301260821748.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/15b9ws1260821748.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/16bwyh1260821748.tab")
+ }
>
> try(system("convert tmp/1xo2f1260821748.ps tmp/1xo2f1260821748.png",intern=TRUE))
character(0)
> try(system("convert tmp/20y4x1260821748.ps tmp/20y4x1260821748.png",intern=TRUE))
character(0)
> try(system("convert tmp/3jdjl1260821748.ps tmp/3jdjl1260821748.png",intern=TRUE))
character(0)
> try(system("convert tmp/4s5qi1260821748.ps tmp/4s5qi1260821748.png",intern=TRUE))
character(0)
> try(system("convert tmp/5doy51260821748.ps tmp/5doy51260821748.png",intern=TRUE))
character(0)
> try(system("convert tmp/62of21260821748.ps tmp/62of21260821748.png",intern=TRUE))
character(0)
> try(system("convert tmp/7la1e1260821748.ps tmp/7la1e1260821748.png",intern=TRUE))
character(0)
> try(system("convert tmp/8a0ml1260821748.ps tmp/8a0ml1260821748.png",intern=TRUE))
character(0)
> try(system("convert tmp/9wpgc1260821748.ps tmp/9wpgc1260821748.png",intern=TRUE))
character(0)
> try(system("convert tmp/106hd31260821748.ps tmp/106hd31260821748.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.344 1.568 2.916