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(29.837,0,29.571,0,30.167,0,30.524,0,30.996,0,31.033,0,31.198,0,30.937,0,31.649,0,33.115,0,34.106,0,33.926,0,33.382,0,32.851,0,32.948,0,36.112,0,36.113,0,35.210,0,35.193,0,34.383,0,35.349,0,37.058,0,38.076,0,36.630,0,36.045,0,35.638,0,35.114,0,35.465,0,35.254,0,35.299,0,35.916,0,36.683,0,37.288,0,38.536,0,38.977,0,36.407,0,34.955,0,34.951,0,32.680,0,34.791,0,34.178,0,35.213,0,34.871,0,35.299,0,35.443,0,37.108,0,36.419,0,34.471,0,33.868,0,34.385,0,33.643,1,34.627,1,32.919,1,35.500,1,36.110,1,37.086,1,37.711,1,40.427,1,39.884,1,38.512,1,38.767,1),dim=c(2,61),dimnames=list(c('saldo_zichtrek','crisis'),1:61))
> y <- array(NA,dim=c(2,61),dimnames=list(c('saldo_zichtrek','crisis'),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 = 'No 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
saldo_zichtrek crisis M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 29.837 0 1 0 0 0 0 0 0 0 0 0 0
2 29.571 0 0 1 0 0 0 0 0 0 0 0 0
3 30.167 0 0 0 1 0 0 0 0 0 0 0 0
4 30.524 0 0 0 0 1 0 0 0 0 0 0 0
5 30.996 0 0 0 0 0 1 0 0 0 0 0 0
6 31.033 0 0 0 0 0 0 1 0 0 0 0 0
7 31.198 0 0 0 0 0 0 0 1 0 0 0 0
8 30.937 0 0 0 0 0 0 0 0 1 0 0 0
9 31.649 0 0 0 0 0 0 0 0 0 1 0 0
10 33.115 0 0 0 0 0 0 0 0 0 0 1 0
11 34.106 0 0 0 0 0 0 0 0 0 0 0 1
12 33.926 0 0 0 0 0 0 0 0 0 0 0 0
13 33.382 0 1 0 0 0 0 0 0 0 0 0 0
14 32.851 0 0 1 0 0 0 0 0 0 0 0 0
15 32.948 0 0 0 1 0 0 0 0 0 0 0 0
16 36.112 0 0 0 0 1 0 0 0 0 0 0 0
17 36.113 0 0 0 0 0 1 0 0 0 0 0 0
18 35.210 0 0 0 0 0 0 1 0 0 0 0 0
19 35.193 0 0 0 0 0 0 0 1 0 0 0 0
20 34.383 0 0 0 0 0 0 0 0 1 0 0 0
21 35.349 0 0 0 0 0 0 0 0 0 1 0 0
22 37.058 0 0 0 0 0 0 0 0 0 0 1 0
23 38.076 0 0 0 0 0 0 0 0 0 0 0 1
24 36.630 0 0 0 0 0 0 0 0 0 0 0 0
25 36.045 0 1 0 0 0 0 0 0 0 0 0 0
26 35.638 0 0 1 0 0 0 0 0 0 0 0 0
27 35.114 0 0 0 1 0 0 0 0 0 0 0 0
28 35.465 0 0 0 0 1 0 0 0 0 0 0 0
29 35.254 0 0 0 0 0 1 0 0 0 0 0 0
30 35.299 0 0 0 0 0 0 1 0 0 0 0 0
31 35.916 0 0 0 0 0 0 0 1 0 0 0 0
32 36.683 0 0 0 0 0 0 0 0 1 0 0 0
33 37.288 0 0 0 0 0 0 0 0 0 1 0 0
34 38.536 0 0 0 0 0 0 0 0 0 0 1 0
35 38.977 0 0 0 0 0 0 0 0 0 0 0 1
36 36.407 0 0 0 0 0 0 0 0 0 0 0 0
37 34.955 0 1 0 0 0 0 0 0 0 0 0 0
38 34.951 0 0 1 0 0 0 0 0 0 0 0 0
39 32.680 0 0 0 1 0 0 0 0 0 0 0 0
40 34.791 0 0 0 0 1 0 0 0 0 0 0 0
41 34.178 0 0 0 0 0 1 0 0 0 0 0 0
42 35.213 0 0 0 0 0 0 1 0 0 0 0 0
43 34.871 0 0 0 0 0 0 0 1 0 0 0 0
44 35.299 0 0 0 0 0 0 0 0 1 0 0 0
45 35.443 0 0 0 0 0 0 0 0 0 1 0 0
46 37.108 0 0 0 0 0 0 0 0 0 0 1 0
47 36.419 0 0 0 0 0 0 0 0 0 0 0 1
48 34.471 0 0 0 0 0 0 0 0 0 0 0 0
49 33.868 0 1 0 0 0 0 0 0 0 0 0 0
50 34.385 0 0 1 0 0 0 0 0 0 0 0 0
51 33.643 1 0 0 1 0 0 0 0 0 0 0 0
52 34.627 1 0 0 0 1 0 0 0 0 0 0 0
53 32.919 1 0 0 0 0 1 0 0 0 0 0 0
54 35.500 1 0 0 0 0 0 1 0 0 0 0 0
55 36.110 1 0 0 0 0 0 0 1 0 0 0 0
56 37.086 1 0 0 0 0 0 0 0 1 0 0 0
57 37.711 1 0 0 0 0 0 0 0 0 1 0 0
58 40.427 1 0 0 0 0 0 0 0 0 0 1 0
59 39.884 1 0 0 0 0 0 0 0 0 0 0 1
60 38.512 1 0 0 0 0 0 0 0 0 0 0 0
61 38.767 1 1 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) crisis M1 M2 M3 M4
35.5500 2.1962 -1.4403 -2.0708 -3.0788 -1.6854
M5 M6 M7 M8 M9 M10
-2.0972 -1.5382 -1.3316 -1.1116 -0.5012 1.2596
M11
1.5032
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-4.2726 -1.0243 0.4768 1.2872 2.6602
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 35.5500 0.9556 37.202 < 2e-16 ***
crisis 2.1962 0.7109 3.089 0.00333 **
M1 -1.4403 1.2797 -1.126 0.26596
M2 -2.0708 1.3439 -1.541 0.12992
M3 -3.0788 1.3364 -2.304 0.02560 *
M4 -1.6854 1.3364 -1.261 0.21334
M5 -2.0972 1.3364 -1.569 0.12314
M6 -1.5382 1.3364 -1.151 0.25542
M7 -1.3316 1.3364 -0.996 0.32404
M8 -1.1116 1.3364 -0.832 0.40964
M9 -0.5012 1.3364 -0.375 0.70928
M10 1.2596 1.3364 0.943 0.35063
M11 1.5032 1.3364 1.125 0.26625
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.113 on 48 degrees of freedom
Multiple R-squared: 0.4126, Adjusted R-squared: 0.2657
F-statistic: 2.809 on 12 and 48 DF, p-value: 0.005493
> 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.9993851 0.0012298110 6.149055e-04
[2,] 0.9999124 0.0001752877 8.764386e-05
[3,] 0.9999289 0.0001422419 7.112094e-05
[4,] 0.9999267 0.0001466140 7.330701e-05
[5,] 0.9999282 0.0001436457 7.182283e-05
[6,] 0.9999205 0.0001589472 7.947358e-05
[7,] 0.9999220 0.0001560240 7.801201e-05
[8,] 0.9998964 0.0002071941 1.035970e-04
[9,] 0.9998072 0.0003855558 1.927779e-04
[10,] 0.9998237 0.0003525982 1.762991e-04
[11,] 0.9998308 0.0003383912 1.691956e-04
[12,] 0.9999012 0.0001976932 9.884659e-05
[13,] 0.9998406 0.0003187788 1.593894e-04
[14,] 0.9998615 0.0002769589 1.384794e-04
[15,] 0.9997148 0.0005704005 2.852003e-04
[16,] 0.9995489 0.0009022032 4.511016e-04
[17,] 0.9994567 0.0010866253 5.433127e-04
[18,] 0.9993688 0.0012624726 6.312363e-04
[19,] 0.9988580 0.0022840482 1.142024e-03
[20,] 0.9985606 0.0028788327 1.439416e-03
[21,] 0.9969409 0.0061181956 3.059098e-03
[22,] 0.9934033 0.0131934697 6.596735e-03
[23,] 0.9867992 0.0264016530 1.320083e-02
[24,] 0.9741784 0.0516432996 2.582165e-02
[25,] 0.9670064 0.0659871639 3.299358e-02
[26,] 0.9869936 0.0260128240 1.300641e-02
[27,] 0.9909404 0.0181192354 9.059618e-03
[28,] 0.9902004 0.0195992869 9.799643e-03
[29,] 0.9883464 0.0233071811 1.165359e-02
[30,] 0.9874024 0.0251952961 1.259765e-02
> postscript(file="/var/www/html/rcomp/tmp/1fc3y1258736559.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/2xip11258736559.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/3oa7s1258736559.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/444hy1258736559.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/5w29c1258736559.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 6
-4.27263774 -3.90820000 -2.30416528 -3.34056528 -2.45676528 -2.97876528
7 8 9 10 11 12
-3.02036528 -3.50136528 -3.39976528 -3.69456528 -2.94716528 -1.62396528
13 14 15 16 17 18
-0.72763774 -0.62820000 0.47683472 2.24743472 2.66023472 1.19823472
19 20 21 22 23 24
0.97463472 -0.05536528 0.30023472 0.24843472 1.02283472 1.08003472
25 26 27 28 29 30
1.93536226 2.15880000 2.64283472 1.60043472 1.80123472 1.28723472
31 32 33 34 35 36
1.69763472 2.24463472 2.23923472 1.72643472 1.92383472 0.85703472
37 38 39 40 41 42
0.84536226 1.47180000 0.20883472 0.92643472 0.72523472 1.20123472
43 44 45 46 47 48
0.65263472 0.86063472 0.39423472 0.29843472 -0.63416528 -1.07896528
49 50 51 52 53 54
-0.24163774 0.90580000 -1.02433887 -1.43373887 -2.72993887 -0.70793887
55 56 57 58 59 60
-0.30453887 0.45146113 0.46606113 1.42126113 0.63466113 0.76586113
61
2.46118868
> postscript(file="/var/www/html/rcomp/tmp/67gm11258736559.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 -4.27263774 NA
1 -3.90820000 -4.27263774
2 -2.30416528 -3.90820000
3 -3.34056528 -2.30416528
4 -2.45676528 -3.34056528
5 -2.97876528 -2.45676528
6 -3.02036528 -2.97876528
7 -3.50136528 -3.02036528
8 -3.39976528 -3.50136528
9 -3.69456528 -3.39976528
10 -2.94716528 -3.69456528
11 -1.62396528 -2.94716528
12 -0.72763774 -1.62396528
13 -0.62820000 -0.72763774
14 0.47683472 -0.62820000
15 2.24743472 0.47683472
16 2.66023472 2.24743472
17 1.19823472 2.66023472
18 0.97463472 1.19823472
19 -0.05536528 0.97463472
20 0.30023472 -0.05536528
21 0.24843472 0.30023472
22 1.02283472 0.24843472
23 1.08003472 1.02283472
24 1.93536226 1.08003472
25 2.15880000 1.93536226
26 2.64283472 2.15880000
27 1.60043472 2.64283472
28 1.80123472 1.60043472
29 1.28723472 1.80123472
30 1.69763472 1.28723472
31 2.24463472 1.69763472
32 2.23923472 2.24463472
33 1.72643472 2.23923472
34 1.92383472 1.72643472
35 0.85703472 1.92383472
36 0.84536226 0.85703472
37 1.47180000 0.84536226
38 0.20883472 1.47180000
39 0.92643472 0.20883472
40 0.72523472 0.92643472
41 1.20123472 0.72523472
42 0.65263472 1.20123472
43 0.86063472 0.65263472
44 0.39423472 0.86063472
45 0.29843472 0.39423472
46 -0.63416528 0.29843472
47 -1.07896528 -0.63416528
48 -0.24163774 -1.07896528
49 0.90580000 -0.24163774
50 -1.02433887 0.90580000
51 -1.43373887 -1.02433887
52 -2.72993887 -1.43373887
53 -0.70793887 -2.72993887
54 -0.30453887 -0.70793887
55 0.45146113 -0.30453887
56 0.46606113 0.45146113
57 1.42126113 0.46606113
58 0.63466113 1.42126113
59 0.76586113 0.63466113
60 2.46118868 0.76586113
61 NA 2.46118868
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -3.90820000 -4.27263774
[2,] -2.30416528 -3.90820000
[3,] -3.34056528 -2.30416528
[4,] -2.45676528 -3.34056528
[5,] -2.97876528 -2.45676528
[6,] -3.02036528 -2.97876528
[7,] -3.50136528 -3.02036528
[8,] -3.39976528 -3.50136528
[9,] -3.69456528 -3.39976528
[10,] -2.94716528 -3.69456528
[11,] -1.62396528 -2.94716528
[12,] -0.72763774 -1.62396528
[13,] -0.62820000 -0.72763774
[14,] 0.47683472 -0.62820000
[15,] 2.24743472 0.47683472
[16,] 2.66023472 2.24743472
[17,] 1.19823472 2.66023472
[18,] 0.97463472 1.19823472
[19,] -0.05536528 0.97463472
[20,] 0.30023472 -0.05536528
[21,] 0.24843472 0.30023472
[22,] 1.02283472 0.24843472
[23,] 1.08003472 1.02283472
[24,] 1.93536226 1.08003472
[25,] 2.15880000 1.93536226
[26,] 2.64283472 2.15880000
[27,] 1.60043472 2.64283472
[28,] 1.80123472 1.60043472
[29,] 1.28723472 1.80123472
[30,] 1.69763472 1.28723472
[31,] 2.24463472 1.69763472
[32,] 2.23923472 2.24463472
[33,] 1.72643472 2.23923472
[34,] 1.92383472 1.72643472
[35,] 0.85703472 1.92383472
[36,] 0.84536226 0.85703472
[37,] 1.47180000 0.84536226
[38,] 0.20883472 1.47180000
[39,] 0.92643472 0.20883472
[40,] 0.72523472 0.92643472
[41,] 1.20123472 0.72523472
[42,] 0.65263472 1.20123472
[43,] 0.86063472 0.65263472
[44,] 0.39423472 0.86063472
[45,] 0.29843472 0.39423472
[46,] -0.63416528 0.29843472
[47,] -1.07896528 -0.63416528
[48,] -0.24163774 -1.07896528
[49,] 0.90580000 -0.24163774
[50,] -1.02433887 0.90580000
[51,] -1.43373887 -1.02433887
[52,] -2.72993887 -1.43373887
[53,] -0.70793887 -2.72993887
[54,] -0.30453887 -0.70793887
[55,] 0.45146113 -0.30453887
[56,] 0.46606113 0.45146113
[57,] 1.42126113 0.46606113
[58,] 0.63466113 1.42126113
[59,] 0.76586113 0.63466113
[60,] 2.46118868 0.76586113
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -3.90820000 -4.27263774
2 -2.30416528 -3.90820000
3 -3.34056528 -2.30416528
4 -2.45676528 -3.34056528
5 -2.97876528 -2.45676528
6 -3.02036528 -2.97876528
7 -3.50136528 -3.02036528
8 -3.39976528 -3.50136528
9 -3.69456528 -3.39976528
10 -2.94716528 -3.69456528
11 -1.62396528 -2.94716528
12 -0.72763774 -1.62396528
13 -0.62820000 -0.72763774
14 0.47683472 -0.62820000
15 2.24743472 0.47683472
16 2.66023472 2.24743472
17 1.19823472 2.66023472
18 0.97463472 1.19823472
19 -0.05536528 0.97463472
20 0.30023472 -0.05536528
21 0.24843472 0.30023472
22 1.02283472 0.24843472
23 1.08003472 1.02283472
24 1.93536226 1.08003472
25 2.15880000 1.93536226
26 2.64283472 2.15880000
27 1.60043472 2.64283472
28 1.80123472 1.60043472
29 1.28723472 1.80123472
30 1.69763472 1.28723472
31 2.24463472 1.69763472
32 2.23923472 2.24463472
33 1.72643472 2.23923472
34 1.92383472 1.72643472
35 0.85703472 1.92383472
36 0.84536226 0.85703472
37 1.47180000 0.84536226
38 0.20883472 1.47180000
39 0.92643472 0.20883472
40 0.72523472 0.92643472
41 1.20123472 0.72523472
42 0.65263472 1.20123472
43 0.86063472 0.65263472
44 0.39423472 0.86063472
45 0.29843472 0.39423472
46 -0.63416528 0.29843472
47 -1.07896528 -0.63416528
48 -0.24163774 -1.07896528
49 0.90580000 -0.24163774
50 -1.02433887 0.90580000
51 -1.43373887 -1.02433887
52 -2.72993887 -1.43373887
53 -0.70793887 -2.72993887
54 -0.30453887 -0.70793887
55 0.45146113 -0.30453887
56 0.46606113 0.45146113
57 1.42126113 0.46606113
58 0.63466113 1.42126113
59 0.76586113 0.63466113
60 2.46118868 0.76586113
> 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/7vo8z1258736559.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/87okv1258736559.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/9ixql1258736559.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/10qoic1258736559.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/11s5s01258736559.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/12bwqr1258736559.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/13mzgu1258736559.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/14pma71258736559.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/15nunx1258736559.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/16lik21258736559.tab")
+ }
>
> system("convert tmp/1fc3y1258736559.ps tmp/1fc3y1258736559.png")
> system("convert tmp/2xip11258736559.ps tmp/2xip11258736559.png")
> system("convert tmp/3oa7s1258736559.ps tmp/3oa7s1258736559.png")
> system("convert tmp/444hy1258736559.ps tmp/444hy1258736559.png")
> system("convert tmp/5w29c1258736559.ps tmp/5w29c1258736559.png")
> system("convert tmp/67gm11258736559.ps tmp/67gm11258736559.png")
> system("convert tmp/7vo8z1258736559.ps tmp/7vo8z1258736559.png")
> system("convert tmp/87okv1258736559.ps tmp/87okv1258736559.png")
> system("convert tmp/9ixql1258736559.ps tmp/9ixql1258736559.png")
> system("convert tmp/10qoic1258736559.ps tmp/10qoic1258736559.png")
>
>
> proc.time()
user system elapsed
2.377 1.534 2.887