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(562000,4814,561000,3908,555000,5250,544000,3937,537000,4004,543000,5560,594000,3922,611000,3759,613000,4138,611000,4634,594000,3996,595000,4308,591000,4143,589000,4429,584000,5219,573000,4929,567000,5755,569000,5592,621000,4163,629000,4962,628000,5208,612000,4755,595000,4491,597000,5732,593000,5731,590000,5040,580000,6102,574000,4904,573000,5369,573000,5578,620000,4619,626000,4731,620000,5011,588000,5299,566000,4146,557000,4625,561000,4736,549000,4219,532000,5116,526000,4205,511000,4121,499000,5103,555000,4300,565000,4578,542000,3809,527000,5526,510000,4247,514000,3830,517000,4394,508000,4826,493000,4409,490000,4569,469000,4106,478000,4794,528000,3914,534000,3793,518000,4405,506000,4022,502000,4100,516000,4788),dim=c(2,60),dimnames=list(c('werkloos','bouw'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('werkloos','bouw'),1:60))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
werkloos bouw M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 562000 4814 1 0 0 0 0 0 0 0 0 0 0 1
2 561000 3908 0 1 0 0 0 0 0 0 0 0 0 2
3 555000 5250 0 0 1 0 0 0 0 0 0 0 0 3
4 544000 3937 0 0 0 1 0 0 0 0 0 0 0 4
5 537000 4004 0 0 0 0 1 0 0 0 0 0 0 5
6 543000 5560 0 0 0 0 0 1 0 0 0 0 0 6
7 594000 3922 0 0 0 0 0 0 1 0 0 0 0 7
8 611000 3759 0 0 0 0 0 0 0 1 0 0 0 8
9 613000 4138 0 0 0 0 0 0 0 0 1 0 0 9
10 611000 4634 0 0 0 0 0 0 0 0 0 1 0 10
11 594000 3996 0 0 0 0 0 0 0 0 0 0 1 11
12 595000 4308 0 0 0 0 0 0 0 0 0 0 0 12
13 591000 4143 1 0 0 0 0 0 0 0 0 0 0 13
14 589000 4429 0 1 0 0 0 0 0 0 0 0 0 14
15 584000 5219 0 0 1 0 0 0 0 0 0 0 0 15
16 573000 4929 0 0 0 1 0 0 0 0 0 0 0 16
17 567000 5755 0 0 0 0 1 0 0 0 0 0 0 17
18 569000 5592 0 0 0 0 0 1 0 0 0 0 0 18
19 621000 4163 0 0 0 0 0 0 1 0 0 0 0 19
20 629000 4962 0 0 0 0 0 0 0 1 0 0 0 20
21 628000 5208 0 0 0 0 0 0 0 0 1 0 0 21
22 612000 4755 0 0 0 0 0 0 0 0 0 1 0 22
23 595000 4491 0 0 0 0 0 0 0 0 0 0 1 23
24 597000 5732 0 0 0 0 0 0 0 0 0 0 0 24
25 593000 5731 1 0 0 0 0 0 0 0 0 0 0 25
26 590000 5040 0 1 0 0 0 0 0 0 0 0 0 26
27 580000 6102 0 0 1 0 0 0 0 0 0 0 0 27
28 574000 4904 0 0 0 1 0 0 0 0 0 0 0 28
29 573000 5369 0 0 0 0 1 0 0 0 0 0 0 29
30 573000 5578 0 0 0 0 0 1 0 0 0 0 0 30
31 620000 4619 0 0 0 0 0 0 1 0 0 0 0 31
32 626000 4731 0 0 0 0 0 0 0 1 0 0 0 32
33 620000 5011 0 0 0 0 0 0 0 0 1 0 0 33
34 588000 5299 0 0 0 0 0 0 0 0 0 1 0 34
35 566000 4146 0 0 0 0 0 0 0 0 0 0 1 35
36 557000 4625 0 0 0 0 0 0 0 0 0 0 0 36
37 561000 4736 1 0 0 0 0 0 0 0 0 0 0 37
38 549000 4219 0 1 0 0 0 0 0 0 0 0 0 38
39 532000 5116 0 0 1 0 0 0 0 0 0 0 0 39
40 526000 4205 0 0 0 1 0 0 0 0 0 0 0 40
41 511000 4121 0 0 0 0 1 0 0 0 0 0 0 41
42 499000 5103 0 0 0 0 0 1 0 0 0 0 0 42
43 555000 4300 0 0 0 0 0 0 1 0 0 0 0 43
44 565000 4578 0 0 0 0 0 0 0 1 0 0 0 44
45 542000 3809 0 0 0 0 0 0 0 0 1 0 0 45
46 527000 5526 0 0 0 0 0 0 0 0 0 1 0 46
47 510000 4247 0 0 0 0 0 0 0 0 0 0 1 47
48 514000 3830 0 0 0 0 0 0 0 0 0 0 0 48
49 517000 4394 1 0 0 0 0 0 0 0 0 0 0 49
50 508000 4826 0 1 0 0 0 0 0 0 0 0 0 50
51 493000 4409 0 0 1 0 0 0 0 0 0 0 0 51
52 490000 4569 0 0 0 1 0 0 0 0 0 0 0 52
53 469000 4106 0 0 0 0 1 0 0 0 0 0 0 53
54 478000 4794 0 0 0 0 0 1 0 0 0 0 0 54
55 528000 3914 0 0 0 0 0 0 1 0 0 0 0 55
56 534000 3793 0 0 0 0 0 0 0 1 0 0 0 56
57 518000 4405 0 0 0 0 0 0 0 0 1 0 0 57
58 506000 4022 0 0 0 0 0 0 0 0 0 1 0 58
59 502000 4100 0 0 0 0 0 0 0 0 0 0 1 59
60 516000 4788 0 0 0 0 0 0 0 0 0 0 0 60
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) bouw M1 M2 M3 M4
475949.68 29.57 -11842.49 -7378.79 -38101.41 -22886.12
M5 M6 M7 M8 M9 M10
-36075.64 -52820.64 33752.13 39406.66 27789.76 4149.29
M11 t
9613.89 -1607.11
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-42861 -10682 1528 13713 33130
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 475949.679 27487.145 17.315 < 2e-16 ***
bouw 29.572 5.301 5.579 1.23e-06 ***
M1 -11842.490 12629.736 -0.938 0.353313
M2 -7378.787 12656.746 -0.583 0.562746
M3 -38101.407 12894.808 -2.955 0.004919 **
M4 -22886.124 12610.279 -1.815 0.076066 .
M5 -36075.639 12560.794 -2.872 0.006148 **
M6 -52820.645 13005.695 -4.061 0.000188 ***
M7 33752.133 12808.504 2.635 0.011424 *
M8 39406.659 12636.769 3.118 0.003134 **
M9 27789.755 12551.186 2.214 0.031817 *
M10 4149.290 12557.245 0.330 0.742575
M11 9613.892 12757.193 0.754 0.454926
t -1607.114 152.058 -10.569 6.80e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 19790 on 46 degrees of freedom
Multiple R-squared: 0.8267, Adjusted R-squared: 0.7777
F-statistic: 16.87 on 13 and 46 DF, p-value: 2.531e-13
> 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,] 8.727004e-05 1.745401e-04 0.99991273
[2,] 7.820891e-05 1.564178e-04 0.99992179
[3,] 8.691267e-06 1.738253e-05 0.99999131
[4,] 7.465190e-04 1.493038e-03 0.99925348
[5,] 1.835493e-03 3.670986e-03 0.99816451
[6,] 4.144536e-02 8.289073e-02 0.95855464
[7,] 8.191770e-02 1.638354e-01 0.91808230
[8,] 1.090094e-01 2.180187e-01 0.89099065
[9,] 9.975986e-02 1.995197e-01 0.90024014
[10,] 8.038892e-02 1.607778e-01 0.91961108
[11,] 7.510130e-02 1.502026e-01 0.92489870
[12,] 5.130446e-02 1.026089e-01 0.94869554
[13,] 3.165858e-02 6.331717e-02 0.96834142
[14,] 2.931792e-02 5.863585e-02 0.97068208
[15,] 2.861324e-02 5.722647e-02 0.97138676
[16,] 3.871932e-02 7.743864e-02 0.96128068
[17,] 1.750037e-01 3.500073e-01 0.82499634
[18,] 6.214169e-01 7.571662e-01 0.37858309
[19,] 8.282376e-01 3.435249e-01 0.17176244
[20,] 8.637027e-01 2.725947e-01 0.13629734
[21,] 8.646654e-01 2.706693e-01 0.13533463
[22,] 8.910805e-01 2.178391e-01 0.10891955
[23,] 8.777557e-01 2.444885e-01 0.12224427
[24,] 8.771359e-01 2.457282e-01 0.12286409
[25,] 9.418840e-01 1.162321e-01 0.05811604
[26,] 8.967378e-01 2.065244e-01 0.10326219
[27,] 8.421789e-01 3.156423e-01 0.15782115
> postscript(file="/var/www/html/rcomp/tmp/1m2lz1258815402.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/21rwq1258815402.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/3ux7g1258815402.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/4997k1258815403.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/5urhj1258815403.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 = 60
Frequency = 1
1 2 3 4 5 6
-42861.1781 -19925.2538 -33281.5608 -19061.2865 -13246.0025 -34908.3977
7 8 9 10 11 12
-20434.6172 -2661.7424 1354.3696 9934.0830 7943.7285 10938.1738
13 14 15 16 17 18
25267.2091 11952.9392 15920.5476 -111.6518 -15741.7515 9430.6552
19 20 21 22 23 24
18723.8228 -951.8652 3997.3641 26641.2002 13590.8015 -9887.4297
25 26 27 28 29 30
-408.2533 14169.6242 5093.5642 20913.0228 20958.5272 33130.0343
31 32 33 34 35 36
23524.2160 22164.7054 21108.4761 5839.2301 14078.6155 2134.4850
37 38 39 40 41 42
16301.5626 16733.8581 5537.2293 12869.4348 15150.1377 -7537.7513
43 44 45 46 47 48
-12756.8501 -15025.3643 -2060.2397 -42588.3176 -25622.8210 1929.8388
49 50 51 52 53 54
1700.6597 -22931.1676 6730.2197 -14609.5194 -7120.9108 -114.5405
55 56 57 58 59 60
-9056.5714 -3525.7336 -24399.9700 173.8044 -9990.3245 -5115.0679
> postscript(file="/var/www/html/rcomp/tmp/6gxsi1258815403.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -42861.1781 NA
1 -19925.2538 -42861.1781
2 -33281.5608 -19925.2538
3 -19061.2865 -33281.5608
4 -13246.0025 -19061.2865
5 -34908.3977 -13246.0025
6 -20434.6172 -34908.3977
7 -2661.7424 -20434.6172
8 1354.3696 -2661.7424
9 9934.0830 1354.3696
10 7943.7285 9934.0830
11 10938.1738 7943.7285
12 25267.2091 10938.1738
13 11952.9392 25267.2091
14 15920.5476 11952.9392
15 -111.6518 15920.5476
16 -15741.7515 -111.6518
17 9430.6552 -15741.7515
18 18723.8228 9430.6552
19 -951.8652 18723.8228
20 3997.3641 -951.8652
21 26641.2002 3997.3641
22 13590.8015 26641.2002
23 -9887.4297 13590.8015
24 -408.2533 -9887.4297
25 14169.6242 -408.2533
26 5093.5642 14169.6242
27 20913.0228 5093.5642
28 20958.5272 20913.0228
29 33130.0343 20958.5272
30 23524.2160 33130.0343
31 22164.7054 23524.2160
32 21108.4761 22164.7054
33 5839.2301 21108.4761
34 14078.6155 5839.2301
35 2134.4850 14078.6155
36 16301.5626 2134.4850
37 16733.8581 16301.5626
38 5537.2293 16733.8581
39 12869.4348 5537.2293
40 15150.1377 12869.4348
41 -7537.7513 15150.1377
42 -12756.8501 -7537.7513
43 -15025.3643 -12756.8501
44 -2060.2397 -15025.3643
45 -42588.3176 -2060.2397
46 -25622.8210 -42588.3176
47 1929.8388 -25622.8210
48 1700.6597 1929.8388
49 -22931.1676 1700.6597
50 6730.2197 -22931.1676
51 -14609.5194 6730.2197
52 -7120.9108 -14609.5194
53 -114.5405 -7120.9108
54 -9056.5714 -114.5405
55 -3525.7336 -9056.5714
56 -24399.9700 -3525.7336
57 173.8044 -24399.9700
58 -9990.3245 173.8044
59 -5115.0679 -9990.3245
60 NA -5115.0679
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -19925.2538 -42861.1781
[2,] -33281.5608 -19925.2538
[3,] -19061.2865 -33281.5608
[4,] -13246.0025 -19061.2865
[5,] -34908.3977 -13246.0025
[6,] -20434.6172 -34908.3977
[7,] -2661.7424 -20434.6172
[8,] 1354.3696 -2661.7424
[9,] 9934.0830 1354.3696
[10,] 7943.7285 9934.0830
[11,] 10938.1738 7943.7285
[12,] 25267.2091 10938.1738
[13,] 11952.9392 25267.2091
[14,] 15920.5476 11952.9392
[15,] -111.6518 15920.5476
[16,] -15741.7515 -111.6518
[17,] 9430.6552 -15741.7515
[18,] 18723.8228 9430.6552
[19,] -951.8652 18723.8228
[20,] 3997.3641 -951.8652
[21,] 26641.2002 3997.3641
[22,] 13590.8015 26641.2002
[23,] -9887.4297 13590.8015
[24,] -408.2533 -9887.4297
[25,] 14169.6242 -408.2533
[26,] 5093.5642 14169.6242
[27,] 20913.0228 5093.5642
[28,] 20958.5272 20913.0228
[29,] 33130.0343 20958.5272
[30,] 23524.2160 33130.0343
[31,] 22164.7054 23524.2160
[32,] 21108.4761 22164.7054
[33,] 5839.2301 21108.4761
[34,] 14078.6155 5839.2301
[35,] 2134.4850 14078.6155
[36,] 16301.5626 2134.4850
[37,] 16733.8581 16301.5626
[38,] 5537.2293 16733.8581
[39,] 12869.4348 5537.2293
[40,] 15150.1377 12869.4348
[41,] -7537.7513 15150.1377
[42,] -12756.8501 -7537.7513
[43,] -15025.3643 -12756.8501
[44,] -2060.2397 -15025.3643
[45,] -42588.3176 -2060.2397
[46,] -25622.8210 -42588.3176
[47,] 1929.8388 -25622.8210
[48,] 1700.6597 1929.8388
[49,] -22931.1676 1700.6597
[50,] 6730.2197 -22931.1676
[51,] -14609.5194 6730.2197
[52,] -7120.9108 -14609.5194
[53,] -114.5405 -7120.9108
[54,] -9056.5714 -114.5405
[55,] -3525.7336 -9056.5714
[56,] -24399.9700 -3525.7336
[57,] 173.8044 -24399.9700
[58,] -9990.3245 173.8044
[59,] -5115.0679 -9990.3245
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -19925.2538 -42861.1781
2 -33281.5608 -19925.2538
3 -19061.2865 -33281.5608
4 -13246.0025 -19061.2865
5 -34908.3977 -13246.0025
6 -20434.6172 -34908.3977
7 -2661.7424 -20434.6172
8 1354.3696 -2661.7424
9 9934.0830 1354.3696
10 7943.7285 9934.0830
11 10938.1738 7943.7285
12 25267.2091 10938.1738
13 11952.9392 25267.2091
14 15920.5476 11952.9392
15 -111.6518 15920.5476
16 -15741.7515 -111.6518
17 9430.6552 -15741.7515
18 18723.8228 9430.6552
19 -951.8652 18723.8228
20 3997.3641 -951.8652
21 26641.2002 3997.3641
22 13590.8015 26641.2002
23 -9887.4297 13590.8015
24 -408.2533 -9887.4297
25 14169.6242 -408.2533
26 5093.5642 14169.6242
27 20913.0228 5093.5642
28 20958.5272 20913.0228
29 33130.0343 20958.5272
30 23524.2160 33130.0343
31 22164.7054 23524.2160
32 21108.4761 22164.7054
33 5839.2301 21108.4761
34 14078.6155 5839.2301
35 2134.4850 14078.6155
36 16301.5626 2134.4850
37 16733.8581 16301.5626
38 5537.2293 16733.8581
39 12869.4348 5537.2293
40 15150.1377 12869.4348
41 -7537.7513 15150.1377
42 -12756.8501 -7537.7513
43 -15025.3643 -12756.8501
44 -2060.2397 -15025.3643
45 -42588.3176 -2060.2397
46 -25622.8210 -42588.3176
47 1929.8388 -25622.8210
48 1700.6597 1929.8388
49 -22931.1676 1700.6597
50 6730.2197 -22931.1676
51 -14609.5194 6730.2197
52 -7120.9108 -14609.5194
53 -114.5405 -7120.9108
54 -9056.5714 -114.5405
55 -3525.7336 -9056.5714
56 -24399.9700 -3525.7336
57 173.8044 -24399.9700
58 -9990.3245 173.8044
59 -5115.0679 -9990.3245
> 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/755kn1258815403.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/8khp91258815403.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/9cmre1258815403.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/10dv2l1258815403.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/11k63s1258815403.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/1287c41258815403.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/13wkj81258815403.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/14mg6x1258815403.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/15dtjy1258815403.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/1640f31258815403.tab")
+ }
>
> system("convert tmp/1m2lz1258815402.ps tmp/1m2lz1258815402.png")
> system("convert tmp/21rwq1258815402.ps tmp/21rwq1258815402.png")
> system("convert tmp/3ux7g1258815402.ps tmp/3ux7g1258815402.png")
> system("convert tmp/4997k1258815403.ps tmp/4997k1258815403.png")
> system("convert tmp/5urhj1258815403.ps tmp/5urhj1258815403.png")
> system("convert tmp/6gxsi1258815403.ps tmp/6gxsi1258815403.png")
> system("convert tmp/755kn1258815403.ps tmp/755kn1258815403.png")
> system("convert tmp/8khp91258815403.ps tmp/8khp91258815403.png")
> system("convert tmp/9cmre1258815403.ps tmp/9cmre1258815403.png")
> system("convert tmp/10dv2l1258815403.ps tmp/10dv2l1258815403.png")
>
>
> proc.time()
user system elapsed
2.408 1.544 2.871