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(16198.9
+ ,16896.2
+ ,0
+ ,16554.2
+ ,16698
+ ,0
+ ,19554.2
+ ,19691.6
+ ,0
+ ,15903.8
+ ,15930.7
+ ,0
+ ,18003.8
+ ,17444.6
+ ,0
+ ,18329.6
+ ,17699.4
+ ,0
+ ,16260.7
+ ,15189.8
+ ,0
+ ,14851.9
+ ,15672.7
+ ,0
+ ,18174.1
+ ,17180.8
+ ,0
+ ,18406.6
+ ,17664.9
+ ,0
+ ,18466.5
+ ,17862.9
+ ,0
+ ,16016.5
+ ,16162.3
+ ,0
+ ,17428.5
+ ,17463.6
+ ,0
+ ,17167.2
+ ,16772.1
+ ,0
+ ,19630
+ ,19106.9
+ ,0
+ ,17183.6
+ ,16721.3
+ ,0
+ ,18344.7
+ ,18161.3
+ ,0
+ ,19301.4
+ ,18509.9
+ ,0
+ ,18147.5
+ ,17802.7
+ ,0
+ ,16192.9
+ ,16409.9
+ ,0
+ ,18374.4
+ ,17967.7
+ ,0
+ ,20515.2
+ ,20286.6
+ ,0
+ ,18957.2
+ ,19537.3
+ ,0
+ ,16471.5
+ ,18021.9
+ ,0
+ ,18746.8
+ ,20194.3
+ ,0
+ ,19009.5
+ ,19049.6
+ ,0
+ ,19211.2
+ ,20244.7
+ ,0
+ ,20547.7
+ ,21473.3
+ ,0
+ ,19325.8
+ ,19673.6
+ ,0
+ ,20605.5
+ ,21053.2
+ ,0
+ ,20056.9
+ ,20159.5
+ ,0
+ ,16141.4
+ ,18203.6
+ ,0
+ ,20359.8
+ ,21289.5
+ ,0
+ ,19711.6
+ ,20432.3
+ ,1
+ ,15638.6
+ ,17180.4
+ ,1
+ ,14384.5
+ ,15816.8
+ ,1
+ ,13855.6
+ ,15071.8
+ ,1
+ ,14308.3
+ ,14521.1
+ ,1
+ ,15290.6
+ ,15668.8
+ ,1
+ ,14423.8
+ ,14346.9
+ ,1
+ ,13779.7
+ ,13881
+ ,1
+ ,15686.3
+ ,15465.9
+ ,1
+ ,14733.8
+ ,14238.2
+ ,1
+ ,12522.5
+ ,13557.7
+ ,1
+ ,16189.4
+ ,16127.6
+ ,1
+ ,16059.1
+ ,16793.9
+ ,1
+ ,16007.1
+ ,16014
+ ,1
+ ,15806.8
+ ,16867.9
+ ,1
+ ,15160
+ ,16014.6
+ ,0
+ ,15692.1
+ ,15878.6
+ ,0
+ ,18908.9
+ ,18664.9
+ ,0
+ ,16969.9
+ ,17962.5
+ ,0
+ ,16997.5
+ ,17332.7
+ ,0
+ ,19858.9
+ ,19542.1
+ ,0
+ ,17681.2
+ ,17203.6
+ ,0)
+ ,dim=c(3
+ ,55)
+ ,dimnames=list(c('uitvoer'
+ ,'invoer'
+ ,'crisis')
+ ,1:55))
> y <- array(NA,dim=c(3,55),dimnames=list(c('uitvoer','invoer','crisis'),1:55))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = '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
uitvoer invoer crisis M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 16198.9 16896.2 0 1 0 0 0 0 0 0 0 0 0 0
2 16554.2 16698.0 0 0 1 0 0 0 0 0 0 0 0 0
3 19554.2 19691.6 0 0 0 1 0 0 0 0 0 0 0 0
4 15903.8 15930.7 0 0 0 0 1 0 0 0 0 0 0 0
5 18003.8 17444.6 0 0 0 0 0 1 0 0 0 0 0 0
6 18329.6 17699.4 0 0 0 0 0 0 1 0 0 0 0 0
7 16260.7 15189.8 0 0 0 0 0 0 0 1 0 0 0 0
8 14851.9 15672.7 0 0 0 0 0 0 0 0 1 0 0 0
9 18174.1 17180.8 0 0 0 0 0 0 0 0 0 1 0 0
10 18406.6 17664.9 0 0 0 0 0 0 0 0 0 0 1 0
11 18466.5 17862.9 0 0 0 0 0 0 0 0 0 0 0 1
12 16016.5 16162.3 0 0 0 0 0 0 0 0 0 0 0 0
13 17428.5 17463.6 0 1 0 0 0 0 0 0 0 0 0 0
14 17167.2 16772.1 0 0 1 0 0 0 0 0 0 0 0 0
15 19630.0 19106.9 0 0 0 1 0 0 0 0 0 0 0 0
16 17183.6 16721.3 0 0 0 0 1 0 0 0 0 0 0 0
17 18344.7 18161.3 0 0 0 0 0 1 0 0 0 0 0 0
18 19301.4 18509.9 0 0 0 0 0 0 1 0 0 0 0 0
19 18147.5 17802.7 0 0 0 0 0 0 0 1 0 0 0 0
20 16192.9 16409.9 0 0 0 0 0 0 0 0 1 0 0 0
21 18374.4 17967.7 0 0 0 0 0 0 0 0 0 1 0 0
22 20515.2 20286.6 0 0 0 0 0 0 0 0 0 0 1 0
23 18957.2 19537.3 0 0 0 0 0 0 0 0 0 0 0 1
24 16471.5 18021.9 0 0 0 0 0 0 0 0 0 0 0 0
25 18746.8 20194.3 0 1 0 0 0 0 0 0 0 0 0 0
26 19009.5 19049.6 0 0 1 0 0 0 0 0 0 0 0 0
27 19211.2 20244.7 0 0 0 1 0 0 0 0 0 0 0 0
28 20547.7 21473.3 0 0 0 0 1 0 0 0 0 0 0 0
29 19325.8 19673.6 0 0 0 0 0 1 0 0 0 0 0 0
30 20605.5 21053.2 0 0 0 0 0 0 1 0 0 0 0 0
31 20056.9 20159.5 0 0 0 0 0 0 0 1 0 0 0 0
32 16141.4 18203.6 0 0 0 0 0 0 0 0 1 0 0 0
33 20359.8 21289.5 0 0 0 0 0 0 0 0 0 1 0 0
34 19711.6 20432.3 1 0 0 0 0 0 0 0 0 0 1 0
35 15638.6 17180.4 1 0 0 0 0 0 0 0 0 0 0 1
36 14384.5 15816.8 1 0 0 0 0 0 0 0 0 0 0 0
37 13855.6 15071.8 1 1 0 0 0 0 0 0 0 0 0 0
38 14308.3 14521.1 1 0 1 0 0 0 0 0 0 0 0 0
39 15290.6 15668.8 1 0 0 1 0 0 0 0 0 0 0 0
40 14423.8 14346.9 1 0 0 0 1 0 0 0 0 0 0 0
41 13779.7 13881.0 1 0 0 0 0 1 0 0 0 0 0 0
42 15686.3 15465.9 1 0 0 0 0 0 1 0 0 0 0 0
43 14733.8 14238.2 1 0 0 0 0 0 0 1 0 0 0 0
44 12522.5 13557.7 1 0 0 0 0 0 0 0 1 0 0 0
45 16189.4 16127.6 1 0 0 0 0 0 0 0 0 1 0 0
46 16059.1 16793.9 1 0 0 0 0 0 0 0 0 0 1 0
47 16007.1 16014.0 1 0 0 0 0 0 0 0 0 0 0 1
48 15806.8 16867.9 1 0 0 0 0 0 0 0 0 0 0 0
49 15160.0 16014.6 0 1 0 0 0 0 0 0 0 0 0 0
50 15692.1 15878.6 0 0 1 0 0 0 0 0 0 0 0 0
51 18908.9 18664.9 0 0 0 1 0 0 0 0 0 0 0 0
52 16969.9 17962.5 0 0 0 0 1 0 0 0 0 0 0 0
53 16997.5 17332.7 0 0 0 0 0 1 0 0 0 0 0 0
54 19858.9 19542.1 0 0 0 0 0 0 1 0 0 0 0 0
55 17681.2 17203.6 0 0 0 0 0 0 0 1 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) invoer crisis M1 M2 M3
3885.4464 0.7349 -1001.2832 5.8095 674.0416 1109.7768
M4 M5 M6 M7 M8 M9
616.8825 892.8245 1509.7493 1257.7076 -437.2239 1307.6929
M10 M11
1476.8239 913.0469
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-784.00 -208.30 62.89 273.81 703.73
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.885e+03 8.298e+02 4.682 3.10e-05 ***
invoer 7.349e-01 4.418e-02 16.633 < 2e-16 ***
crisis -1.001e+03 1.813e+02 -5.523 2.06e-06 ***
M1 5.810e+00 2.983e+02 0.019 0.984556
M2 6.740e+02 3.006e+02 2.242 0.030427 *
M3 1.110e+03 3.022e+02 3.672 0.000688 ***
M4 6.169e+02 2.980e+02 2.070 0.044771 *
M5 8.928e+02 2.980e+02 2.996 0.004621 **
M6 1.510e+03 3.007e+02 5.021 1.05e-05 ***
M7 1.258e+03 2.990e+02 4.207 0.000137 ***
M8 -4.372e+02 3.190e+02 -1.371 0.177957
M9 1.308e+03 3.146e+02 4.157 0.000160 ***
M10 1.477e+03 3.241e+02 4.556 4.61e-05 ***
M11 9.130e+02 3.136e+02 2.912 0.005790 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 439.6 on 41 degrees of freedom
Multiple R-squared: 0.964, Adjusted R-squared: 0.9526
F-statistic: 84.46 on 13 and 41 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.6818467 0.6363067 0.3181533
[2,] 0.5367507 0.9264985 0.4632493
[3,] 0.4644918 0.9289836 0.5355082
[4,] 0.5807073 0.8385853 0.4192927
[5,] 0.5327707 0.9344586 0.4672293
[6,] 0.5093270 0.9813461 0.4906730
[7,] 0.5709251 0.8581498 0.4290749
[8,] 0.6090665 0.7818671 0.3909335
[9,] 0.5004616 0.9990768 0.4995384
[10,] 0.5221428 0.9557144 0.4778572
[11,] 0.6538969 0.6922063 0.3461031
[12,] 0.5785681 0.8428638 0.4214319
[13,] 0.5134401 0.9731198 0.4865599
[14,] 0.4259781 0.8519562 0.5740219
[15,] 0.3270828 0.6541656 0.6729172
[16,] 0.3717443 0.7434887 0.6282557
[17,] 0.3366310 0.6732619 0.6633690
[18,] 0.2975134 0.5950268 0.7024866
[19,] 0.6505312 0.6989376 0.3494688
[20,] 0.6093463 0.7813075 0.3906537
[21,] 0.4572341 0.9144682 0.5427659
[22,] 0.3288392 0.6576783 0.6711608
> postscript(file="/var/www/html/freestat/rcomp/tmp/1c99u1290767420.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/freestat/rcomp/tmp/2c99u1290767420.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/freestat/rcomp/tmp/3mi9x1290767420.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/freestat/rcomp/tmp/4mi9x1290767420.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/freestat/rcomp/tmp/5mi9x1290767420.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 = 55
Frequency = 1
1 2 3 4 5 6 7
-108.89982 -276.18022 88.17177 -305.55386 405.98145 -72.38875 -45.01238
8 9 10 11 12 13 14
-113.75057 355.27223 62.88969 541.06210 253.83225 703.73382 282.36577
15 16 17 18 19 20 21
593.65141 393.25635 220.19871 303.79751 -78.35937 685.50181 -22.69853
22 23 24 25 26 27 28
244.87583 -198.70754 -657.73568 15.31891 450.99483 -661.28592 265.24471
29 30 31 32 33 34 35
89.95182 -261.10239 99.09435 -684.13806 -478.39626 335.48814 -784.00462
36 37 38 39 40 41 42
-122.98632 -110.21623 78.94574 -217.90208 379.61954 -198.04561 -73.06903
43 44 45 46 47 48 49
128.67495 112.38682 145.82256 -643.25366 441.65006 526.88974 -499.93668
50 51 52 53 54 55
-536.12612 197.36483 -732.56674 -518.08637 102.76266 -104.39757
> postscript(file="/var/www/html/freestat/rcomp/tmp/6faqi1290767420.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 = 55
Frequency = 1
lag(myerror, k = 1) myerror
0 -108.89982 NA
1 -276.18022 -108.89982
2 88.17177 -276.18022
3 -305.55386 88.17177
4 405.98145 -305.55386
5 -72.38875 405.98145
6 -45.01238 -72.38875
7 -113.75057 -45.01238
8 355.27223 -113.75057
9 62.88969 355.27223
10 541.06210 62.88969
11 253.83225 541.06210
12 703.73382 253.83225
13 282.36577 703.73382
14 593.65141 282.36577
15 393.25635 593.65141
16 220.19871 393.25635
17 303.79751 220.19871
18 -78.35937 303.79751
19 685.50181 -78.35937
20 -22.69853 685.50181
21 244.87583 -22.69853
22 -198.70754 244.87583
23 -657.73568 -198.70754
24 15.31891 -657.73568
25 450.99483 15.31891
26 -661.28592 450.99483
27 265.24471 -661.28592
28 89.95182 265.24471
29 -261.10239 89.95182
30 99.09435 -261.10239
31 -684.13806 99.09435
32 -478.39626 -684.13806
33 335.48814 -478.39626
34 -784.00462 335.48814
35 -122.98632 -784.00462
36 -110.21623 -122.98632
37 78.94574 -110.21623
38 -217.90208 78.94574
39 379.61954 -217.90208
40 -198.04561 379.61954
41 -73.06903 -198.04561
42 128.67495 -73.06903
43 112.38682 128.67495
44 145.82256 112.38682
45 -643.25366 145.82256
46 441.65006 -643.25366
47 526.88974 441.65006
48 -499.93668 526.88974
49 -536.12612 -499.93668
50 197.36483 -536.12612
51 -732.56674 197.36483
52 -518.08637 -732.56674
53 102.76266 -518.08637
54 -104.39757 102.76266
55 NA -104.39757
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -276.18022 -108.89982
[2,] 88.17177 -276.18022
[3,] -305.55386 88.17177
[4,] 405.98145 -305.55386
[5,] -72.38875 405.98145
[6,] -45.01238 -72.38875
[7,] -113.75057 -45.01238
[8,] 355.27223 -113.75057
[9,] 62.88969 355.27223
[10,] 541.06210 62.88969
[11,] 253.83225 541.06210
[12,] 703.73382 253.83225
[13,] 282.36577 703.73382
[14,] 593.65141 282.36577
[15,] 393.25635 593.65141
[16,] 220.19871 393.25635
[17,] 303.79751 220.19871
[18,] -78.35937 303.79751
[19,] 685.50181 -78.35937
[20,] -22.69853 685.50181
[21,] 244.87583 -22.69853
[22,] -198.70754 244.87583
[23,] -657.73568 -198.70754
[24,] 15.31891 -657.73568
[25,] 450.99483 15.31891
[26,] -661.28592 450.99483
[27,] 265.24471 -661.28592
[28,] 89.95182 265.24471
[29,] -261.10239 89.95182
[30,] 99.09435 -261.10239
[31,] -684.13806 99.09435
[32,] -478.39626 -684.13806
[33,] 335.48814 -478.39626
[34,] -784.00462 335.48814
[35,] -122.98632 -784.00462
[36,] -110.21623 -122.98632
[37,] 78.94574 -110.21623
[38,] -217.90208 78.94574
[39,] 379.61954 -217.90208
[40,] -198.04561 379.61954
[41,] -73.06903 -198.04561
[42,] 128.67495 -73.06903
[43,] 112.38682 128.67495
[44,] 145.82256 112.38682
[45,] -643.25366 145.82256
[46,] 441.65006 -643.25366
[47,] 526.88974 441.65006
[48,] -499.93668 526.88974
[49,] -536.12612 -499.93668
[50,] 197.36483 -536.12612
[51,] -732.56674 197.36483
[52,] -518.08637 -732.56674
[53,] 102.76266 -518.08637
[54,] -104.39757 102.76266
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -276.18022 -108.89982
2 88.17177 -276.18022
3 -305.55386 88.17177
4 405.98145 -305.55386
5 -72.38875 405.98145
6 -45.01238 -72.38875
7 -113.75057 -45.01238
8 355.27223 -113.75057
9 62.88969 355.27223
10 541.06210 62.88969
11 253.83225 541.06210
12 703.73382 253.83225
13 282.36577 703.73382
14 593.65141 282.36577
15 393.25635 593.65141
16 220.19871 393.25635
17 303.79751 220.19871
18 -78.35937 303.79751
19 685.50181 -78.35937
20 -22.69853 685.50181
21 244.87583 -22.69853
22 -198.70754 244.87583
23 -657.73568 -198.70754
24 15.31891 -657.73568
25 450.99483 15.31891
26 -661.28592 450.99483
27 265.24471 -661.28592
28 89.95182 265.24471
29 -261.10239 89.95182
30 99.09435 -261.10239
31 -684.13806 99.09435
32 -478.39626 -684.13806
33 335.48814 -478.39626
34 -784.00462 335.48814
35 -122.98632 -784.00462
36 -110.21623 -122.98632
37 78.94574 -110.21623
38 -217.90208 78.94574
39 379.61954 -217.90208
40 -198.04561 379.61954
41 -73.06903 -198.04561
42 128.67495 -73.06903
43 112.38682 128.67495
44 145.82256 112.38682
45 -643.25366 145.82256
46 441.65006 -643.25366
47 526.88974 441.65006
48 -499.93668 526.88974
49 -536.12612 -499.93668
50 197.36483 -536.12612
51 -732.56674 197.36483
52 -518.08637 -732.56674
53 102.76266 -518.08637
54 -104.39757 102.76266
> 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/78jp31290767420.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/freestat/rcomp/tmp/88jp31290767420.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/freestat/rcomp/tmp/98jp31290767420.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/freestat/rcomp/tmp/10jso51290767420.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/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/114tnu1290767420.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/12pb3h1290767420.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/13llj81290767420.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/14pm0e1290767420.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/15a4y21290767420.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/16e5f81290767420.tab")
+ }
>
> try(system("convert tmp/1c99u1290767420.ps tmp/1c99u1290767420.png",intern=TRUE))
character(0)
> try(system("convert tmp/2c99u1290767420.ps tmp/2c99u1290767420.png",intern=TRUE))
character(0)
> try(system("convert tmp/3mi9x1290767420.ps tmp/3mi9x1290767420.png",intern=TRUE))
character(0)
> try(system("convert tmp/4mi9x1290767420.ps tmp/4mi9x1290767420.png",intern=TRUE))
character(0)
> try(system("convert tmp/5mi9x1290767420.ps tmp/5mi9x1290767420.png",intern=TRUE))
character(0)
> try(system("convert tmp/6faqi1290767420.ps tmp/6faqi1290767420.png",intern=TRUE))
character(0)
> try(system("convert tmp/78jp31290767420.ps tmp/78jp31290767420.png",intern=TRUE))
character(0)
> try(system("convert tmp/88jp31290767420.ps tmp/88jp31290767420.png",intern=TRUE))
character(0)
> try(system("convert tmp/98jp31290767420.ps tmp/98jp31290767420.png",intern=TRUE))
character(0)
> try(system("convert tmp/10jso51290767420.ps tmp/10jso51290767420.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.854 2.511 4.325