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(376.974,0,377.632,0,378.205,0,370.861,0,369.167,0,371.551,0,382.842,0,381.903,0,384.502,0,392.058,0,384.359,0,388.884,0,386.586,0,387.495,0,385.705,0,378.67,0,377.367,0,376.911,0,389.827,0,387.82,0,387.267,0,380.575,0,372.402,0,376.74,0,377.795,0,376.126,0,370.804,0,367.98,0,367.866,0,366.121,0,379.421,0,378.519,0,372.423,0,355.072,0,344.693,0,342.892,0,344.178,0,337.606,0,327.103,0,323.953,0,316.532,0,306.307,0,327.225,0,329.573,0,313.761,0,307.836,0,300.074,0,304.198,0,306.122,0,300.414,0,292.133,0,290.616,0,280.244,1,285.179,1,305.486,1,305.957,1,293.886,1,289.441,1,288.776,1,299.149,1,306.532,1,309.914,1,313.468,1,314.901,1,309.16,1,316.15,1,336.544,1,339.196,1,326.738,1,320.838,1,318.62,1,331.533,1,335.378,1),dim=c(2,73),dimnames=list(c('Maandelijksewerkloosheid','x'),1:73))
> y <- array(NA,dim=c(2,73),dimnames=list(c('Maandelijksewerkloosheid','x'),1:73))
> 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
Maandelijksewerkloosheid x M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 376.974 0 1 0 0 0 0 0 0 0 0 0 0 1
2 377.632 0 0 1 0 0 0 0 0 0 0 0 0 2
3 378.205 0 0 0 1 0 0 0 0 0 0 0 0 3
4 370.861 0 0 0 0 1 0 0 0 0 0 0 0 4
5 369.167 0 0 0 0 0 1 0 0 0 0 0 0 5
6 371.551 0 0 0 0 0 0 1 0 0 0 0 0 6
7 382.842 0 0 0 0 0 0 0 1 0 0 0 0 7
8 381.903 0 0 0 0 0 0 0 0 1 0 0 0 8
9 384.502 0 0 0 0 0 0 0 0 0 1 0 0 9
10 392.058 0 0 0 0 0 0 0 0 0 0 1 0 10
11 384.359 0 0 0 0 0 0 0 0 0 0 0 1 11
12 388.884 0 0 0 0 0 0 0 0 0 0 0 0 12
13 386.586 0 1 0 0 0 0 0 0 0 0 0 0 13
14 387.495 0 0 1 0 0 0 0 0 0 0 0 0 14
15 385.705 0 0 0 1 0 0 0 0 0 0 0 0 15
16 378.670 0 0 0 0 1 0 0 0 0 0 0 0 16
17 377.367 0 0 0 0 0 1 0 0 0 0 0 0 17
18 376.911 0 0 0 0 0 0 1 0 0 0 0 0 18
19 389.827 0 0 0 0 0 0 0 1 0 0 0 0 19
20 387.820 0 0 0 0 0 0 0 0 1 0 0 0 20
21 387.267 0 0 0 0 0 0 0 0 0 1 0 0 21
22 380.575 0 0 0 0 0 0 0 0 0 0 1 0 22
23 372.402 0 0 0 0 0 0 0 0 0 0 0 1 23
24 376.740 0 0 0 0 0 0 0 0 0 0 0 0 24
25 377.795 0 1 0 0 0 0 0 0 0 0 0 0 25
26 376.126 0 0 1 0 0 0 0 0 0 0 0 0 26
27 370.804 0 0 0 1 0 0 0 0 0 0 0 0 27
28 367.980 0 0 0 0 1 0 0 0 0 0 0 0 28
29 367.866 0 0 0 0 0 1 0 0 0 0 0 0 29
30 366.121 0 0 0 0 0 0 1 0 0 0 0 0 30
31 379.421 0 0 0 0 0 0 0 1 0 0 0 0 31
32 378.519 0 0 0 0 0 0 0 0 1 0 0 0 32
33 372.423 0 0 0 0 0 0 0 0 0 1 0 0 33
34 355.072 0 0 0 0 0 0 0 0 0 0 1 0 34
35 344.693 0 0 0 0 0 0 0 0 0 0 0 1 35
36 342.892 0 0 0 0 0 0 0 0 0 0 0 0 36
37 344.178 0 1 0 0 0 0 0 0 0 0 0 0 37
38 337.606 0 0 1 0 0 0 0 0 0 0 0 0 38
39 327.103 0 0 0 1 0 0 0 0 0 0 0 0 39
40 323.953 0 0 0 0 1 0 0 0 0 0 0 0 40
41 316.532 0 0 0 0 0 1 0 0 0 0 0 0 41
42 306.307 0 0 0 0 0 0 1 0 0 0 0 0 42
43 327.225 0 0 0 0 0 0 0 1 0 0 0 0 43
44 329.573 0 0 0 0 0 0 0 0 1 0 0 0 44
45 313.761 0 0 0 0 0 0 0 0 0 1 0 0 45
46 307.836 0 0 0 0 0 0 0 0 0 0 1 0 46
47 300.074 0 0 0 0 0 0 0 0 0 0 0 1 47
48 304.198 0 0 0 0 0 0 0 0 0 0 0 0 48
49 306.122 0 1 0 0 0 0 0 0 0 0 0 0 49
50 300.414 0 0 1 0 0 0 0 0 0 0 0 0 50
51 292.133 0 0 0 1 0 0 0 0 0 0 0 0 51
52 290.616 0 0 0 0 1 0 0 0 0 0 0 0 52
53 280.244 1 0 0 0 0 1 0 0 0 0 0 0 53
54 285.179 1 0 0 0 0 0 1 0 0 0 0 0 54
55 305.486 1 0 0 0 0 0 0 1 0 0 0 0 55
56 305.957 1 0 0 0 0 0 0 0 1 0 0 0 56
57 293.886 1 0 0 0 0 0 0 0 0 1 0 0 57
58 289.441 1 0 0 0 0 0 0 0 0 0 1 0 58
59 288.776 1 0 0 0 0 0 0 0 0 0 0 1 59
60 299.149 1 0 0 0 0 0 0 0 0 0 0 0 60
61 306.532 1 1 0 0 0 0 0 0 0 0 0 0 61
62 309.914 1 0 1 0 0 0 0 0 0 0 0 0 62
63 313.468 1 0 0 1 0 0 0 0 0 0 0 0 63
64 314.901 1 0 0 0 1 0 0 0 0 0 0 0 64
65 309.160 1 0 0 0 0 1 0 0 0 0 0 0 65
66 316.150 1 0 0 0 0 0 1 0 0 0 0 0 66
67 336.544 1 0 0 0 0 0 0 1 0 0 0 0 67
68 339.196 1 0 0 0 0 0 0 0 1 0 0 0 68
69 326.738 1 0 0 0 0 0 0 0 0 1 0 0 69
70 320.838 1 0 0 0 0 0 0 0 0 0 1 0 70
71 318.620 1 0 0 0 0 0 0 0 0 0 0 1 71
72 331.533 1 0 0 0 0 0 0 0 0 0 0 0 72
73 335.378 1 1 0 0 0 0 0 0 0 0 0 0 73
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) x M1 M2 M3 M4
402.6813 9.6305 -0.2321 -6.3168 -8.3896 -10.2404
M5 M6 M7 M8 M9 M10
-14.7309 -12.8617 5.2147 7.0405 1.1974 -2.7067
M11 t
-7.3007 -1.5554
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-34.90 -15.62 3.75 15.09 36.84
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 402.6813 10.0666 40.002 < 2e-16 ***
x 9.6305 8.5482 1.127 0.264
M1 -0.2321 11.2874 -0.021 0.984
M2 -6.3168 11.7493 -0.538 0.593
M3 -8.3896 11.7395 -0.715 0.478
M4 -10.2404 11.7327 -0.873 0.386
M5 -14.7309 11.7643 -1.252 0.215
M6 -12.8617 11.7457 -1.095 0.278
M7 5.2147 11.7300 0.445 0.658
M8 7.0405 11.7171 0.601 0.550
M9 1.1974 11.7071 0.102 0.919
M10 -2.7067 11.6999 -0.231 0.818
M11 -7.3007 11.6956 -0.624 0.535
t -1.5554 0.1833 -8.487 8.26e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 20.25 on 59 degrees of freedom
Multiple R-squared: 0.7281, Adjusted R-squared: 0.6682
F-statistic: 12.15 on 13 and 59 DF, p-value: 2.873e-12
> 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,] 9.002089e-05 1.800418e-04 9.999100e-01
[2,] 2.413103e-05 4.826206e-05 9.999759e-01
[3,] 1.404711e-06 2.809422e-06 9.999986e-01
[4,] 1.229813e-07 2.459626e-07 9.999999e-01
[5,] 9.129012e-08 1.825802e-07 9.999999e-01
[6,] 3.427597e-05 6.855194e-05 9.999657e-01
[7,] 7.324716e-05 1.464943e-04 9.999268e-01
[8,] 7.384441e-05 1.476888e-04 9.999262e-01
[9,] 2.674494e-05 5.348987e-05 9.999733e-01
[10,] 1.198336e-05 2.396672e-05 9.999880e-01
[11,] 8.581917e-06 1.716383e-05 9.999914e-01
[12,] 3.398022e-06 6.796044e-06 9.999966e-01
[13,] 1.608575e-06 3.217149e-06 9.999984e-01
[14,] 9.954984e-07 1.990997e-06 9.999990e-01
[15,] 5.338791e-07 1.067758e-06 9.999995e-01
[16,] 3.065088e-07 6.130177e-07 9.999997e-01
[17,] 7.227128e-07 1.445426e-06 9.999993e-01
[18,] 4.340638e-05 8.681276e-05 9.999566e-01
[19,] 7.487390e-04 1.497478e-03 9.992513e-01
[20,] 7.949688e-03 1.589938e-02 9.920503e-01
[21,] 2.835185e-02 5.670369e-02 9.716482e-01
[22,] 1.327260e-01 2.654521e-01 8.672740e-01
[23,] 4.326117e-01 8.652234e-01 5.673883e-01
[24,] 8.094545e-01 3.810910e-01 1.905455e-01
[25,] 9.428233e-01 1.143534e-01 5.717672e-02
[26,] 9.725160e-01 5.496791e-02 2.748396e-02
[27,] 9.818976e-01 3.620476e-02 1.810238e-02
[28,] 9.900301e-01 1.993971e-02 9.969855e-03
[29,] 9.963329e-01 7.334208e-03 3.667104e-03
[30,] 9.993068e-01 1.386358e-03 6.931790e-04
[31,] 9.998780e-01 2.440104e-04 1.220052e-04
[32,] 9.999702e-01 5.958946e-05 2.979473e-05
[33,] 9.999996e-01 7.347830e-07 3.673915e-07
[34,] 1.000000e+00 1.901373e-09 9.506864e-10
[35,] 1.000000e+00 1.160894e-08 5.804472e-09
[36,] 9.999999e-01 1.681771e-07 8.408855e-08
[37,] 9.999995e-01 1.043386e-06 5.216931e-07
[38,] 9.999922e-01 1.560311e-05 7.801554e-06
[39,] 9.998929e-01 2.141355e-04 1.070678e-04
[40,] 9.991326e-01 1.734871e-03 8.674354e-04
> postscript(file="/var/www/html/freestat/rcomp/tmp/163xc1291053163.ps",horizontal=F,onefile=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/2zcwx1291053163.ps",horizontal=F,onefile=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/3zcwx1291053163.ps",horizontal=F,onefile=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/4zcwx1291053163.ps",horizontal=F,onefile=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/5zcwx1291053163.ps",horizontal=F,onefile=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 = 73
Frequency = 1
1 2 3 4 5 6
-23.9198017 -15.6217742 -11.4206075 -15.3584409 -11.0065265 -8.9363599
7 8 9 10 11 12
-14.1663599 -15.3758599 -5.3783599 7.6371401 6.0874735 4.8671401
13 14 15 16 17 18
4.3566071 12.9056345 14.7438012 11.1149679 15.8578822 15.0880489
19 20 21 22 23 24
11.4830489 9.2055489 16.0510489 14.8185489 12.7948822 11.3875489
25 26 27 28 29 30
14.2300158 20.2010433 18.5072100 19.0893766 25.0212910 22.9624576
31 32 33 34 35 36
19.7414576 18.5689576 19.8714576 7.9799576 3.7502910 -3.7960424
37 38 39 40 41 42
-0.7225754 0.3454521 -6.5293813 -6.2732146 -7.6483003 -18.1871336
43 44 45 46 47 48
-13.7901336 -11.7126336 -20.1261336 -20.5916336 -22.2043003 -23.8256336
49 50 51 52 53 54
-20.1141667 -18.1821392 -22.8349725 -20.9458059 -34.9023776 -30.2812109
55 56 57 58 59 60
-26.4952109 -26.2947109 -30.9672109 -29.9527109 -24.4683776 -19.8407109
61 62 63 64 65 66
-10.6702440 0.3517835 7.5339502 12.3731168 12.6780312 19.3541978
67 68 69 70 71 72
23.2271978 25.6086978 20.5491978 20.1086978 24.0400312 31.2076978
73
36.8401648
> postscript(file="/var/www/html/freestat/rcomp/tmp/6ame01291053163.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> dum <- cbind(lag(myerror,k=1),myerror)
> dum
Time Series:
Start = 0
End = 73
Frequency = 1
lag(myerror, k = 1) myerror
0 -23.9198017 NA
1 -15.6217742 -23.9198017
2 -11.4206075 -15.6217742
3 -15.3584409 -11.4206075
4 -11.0065265 -15.3584409
5 -8.9363599 -11.0065265
6 -14.1663599 -8.9363599
7 -15.3758599 -14.1663599
8 -5.3783599 -15.3758599
9 7.6371401 -5.3783599
10 6.0874735 7.6371401
11 4.8671401 6.0874735
12 4.3566071 4.8671401
13 12.9056345 4.3566071
14 14.7438012 12.9056345
15 11.1149679 14.7438012
16 15.8578822 11.1149679
17 15.0880489 15.8578822
18 11.4830489 15.0880489
19 9.2055489 11.4830489
20 16.0510489 9.2055489
21 14.8185489 16.0510489
22 12.7948822 14.8185489
23 11.3875489 12.7948822
24 14.2300158 11.3875489
25 20.2010433 14.2300158
26 18.5072100 20.2010433
27 19.0893766 18.5072100
28 25.0212910 19.0893766
29 22.9624576 25.0212910
30 19.7414576 22.9624576
31 18.5689576 19.7414576
32 19.8714576 18.5689576
33 7.9799576 19.8714576
34 3.7502910 7.9799576
35 -3.7960424 3.7502910
36 -0.7225754 -3.7960424
37 0.3454521 -0.7225754
38 -6.5293813 0.3454521
39 -6.2732146 -6.5293813
40 -7.6483003 -6.2732146
41 -18.1871336 -7.6483003
42 -13.7901336 -18.1871336
43 -11.7126336 -13.7901336
44 -20.1261336 -11.7126336
45 -20.5916336 -20.1261336
46 -22.2043003 -20.5916336
47 -23.8256336 -22.2043003
48 -20.1141667 -23.8256336
49 -18.1821392 -20.1141667
50 -22.8349725 -18.1821392
51 -20.9458059 -22.8349725
52 -34.9023776 -20.9458059
53 -30.2812109 -34.9023776
54 -26.4952109 -30.2812109
55 -26.2947109 -26.4952109
56 -30.9672109 -26.2947109
57 -29.9527109 -30.9672109
58 -24.4683776 -29.9527109
59 -19.8407109 -24.4683776
60 -10.6702440 -19.8407109
61 0.3517835 -10.6702440
62 7.5339502 0.3517835
63 12.3731168 7.5339502
64 12.6780312 12.3731168
65 19.3541978 12.6780312
66 23.2271978 19.3541978
67 25.6086978 23.2271978
68 20.5491978 25.6086978
69 20.1086978 20.5491978
70 24.0400312 20.1086978
71 31.2076978 24.0400312
72 36.8401648 31.2076978
73 NA 36.8401648
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -15.6217742 -23.9198017
[2,] -11.4206075 -15.6217742
[3,] -15.3584409 -11.4206075
[4,] -11.0065265 -15.3584409
[5,] -8.9363599 -11.0065265
[6,] -14.1663599 -8.9363599
[7,] -15.3758599 -14.1663599
[8,] -5.3783599 -15.3758599
[9,] 7.6371401 -5.3783599
[10,] 6.0874735 7.6371401
[11,] 4.8671401 6.0874735
[12,] 4.3566071 4.8671401
[13,] 12.9056345 4.3566071
[14,] 14.7438012 12.9056345
[15,] 11.1149679 14.7438012
[16,] 15.8578822 11.1149679
[17,] 15.0880489 15.8578822
[18,] 11.4830489 15.0880489
[19,] 9.2055489 11.4830489
[20,] 16.0510489 9.2055489
[21,] 14.8185489 16.0510489
[22,] 12.7948822 14.8185489
[23,] 11.3875489 12.7948822
[24,] 14.2300158 11.3875489
[25,] 20.2010433 14.2300158
[26,] 18.5072100 20.2010433
[27,] 19.0893766 18.5072100
[28,] 25.0212910 19.0893766
[29,] 22.9624576 25.0212910
[30,] 19.7414576 22.9624576
[31,] 18.5689576 19.7414576
[32,] 19.8714576 18.5689576
[33,] 7.9799576 19.8714576
[34,] 3.7502910 7.9799576
[35,] -3.7960424 3.7502910
[36,] -0.7225754 -3.7960424
[37,] 0.3454521 -0.7225754
[38,] -6.5293813 0.3454521
[39,] -6.2732146 -6.5293813
[40,] -7.6483003 -6.2732146
[41,] -18.1871336 -7.6483003
[42,] -13.7901336 -18.1871336
[43,] -11.7126336 -13.7901336
[44,] -20.1261336 -11.7126336
[45,] -20.5916336 -20.1261336
[46,] -22.2043003 -20.5916336
[47,] -23.8256336 -22.2043003
[48,] -20.1141667 -23.8256336
[49,] -18.1821392 -20.1141667
[50,] -22.8349725 -18.1821392
[51,] -20.9458059 -22.8349725
[52,] -34.9023776 -20.9458059
[53,] -30.2812109 -34.9023776
[54,] -26.4952109 -30.2812109
[55,] -26.2947109 -26.4952109
[56,] -30.9672109 -26.2947109
[57,] -29.9527109 -30.9672109
[58,] -24.4683776 -29.9527109
[59,] -19.8407109 -24.4683776
[60,] -10.6702440 -19.8407109
[61,] 0.3517835 -10.6702440
[62,] 7.5339502 0.3517835
[63,] 12.3731168 7.5339502
[64,] 12.6780312 12.3731168
[65,] 19.3541978 12.6780312
[66,] 23.2271978 19.3541978
[67,] 25.6086978 23.2271978
[68,] 20.5491978 25.6086978
[69,] 20.1086978 20.5491978
[70,] 24.0400312 20.1086978
[71,] 31.2076978 24.0400312
[72,] 36.8401648 31.2076978
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -15.6217742 -23.9198017
2 -11.4206075 -15.6217742
3 -15.3584409 -11.4206075
4 -11.0065265 -15.3584409
5 -8.9363599 -11.0065265
6 -14.1663599 -8.9363599
7 -15.3758599 -14.1663599
8 -5.3783599 -15.3758599
9 7.6371401 -5.3783599
10 6.0874735 7.6371401
11 4.8671401 6.0874735
12 4.3566071 4.8671401
13 12.9056345 4.3566071
14 14.7438012 12.9056345
15 11.1149679 14.7438012
16 15.8578822 11.1149679
17 15.0880489 15.8578822
18 11.4830489 15.0880489
19 9.2055489 11.4830489
20 16.0510489 9.2055489
21 14.8185489 16.0510489
22 12.7948822 14.8185489
23 11.3875489 12.7948822
24 14.2300158 11.3875489
25 20.2010433 14.2300158
26 18.5072100 20.2010433
27 19.0893766 18.5072100
28 25.0212910 19.0893766
29 22.9624576 25.0212910
30 19.7414576 22.9624576
31 18.5689576 19.7414576
32 19.8714576 18.5689576
33 7.9799576 19.8714576
34 3.7502910 7.9799576
35 -3.7960424 3.7502910
36 -0.7225754 -3.7960424
37 0.3454521 -0.7225754
38 -6.5293813 0.3454521
39 -6.2732146 -6.5293813
40 -7.6483003 -6.2732146
41 -18.1871336 -7.6483003
42 -13.7901336 -18.1871336
43 -11.7126336 -13.7901336
44 -20.1261336 -11.7126336
45 -20.5916336 -20.1261336
46 -22.2043003 -20.5916336
47 -23.8256336 -22.2043003
48 -20.1141667 -23.8256336
49 -18.1821392 -20.1141667
50 -22.8349725 -18.1821392
51 -20.9458059 -22.8349725
52 -34.9023776 -20.9458059
53 -30.2812109 -34.9023776
54 -26.4952109 -30.2812109
55 -26.2947109 -26.4952109
56 -30.9672109 -26.2947109
57 -29.9527109 -30.9672109
58 -24.4683776 -29.9527109
59 -19.8407109 -24.4683776
60 -10.6702440 -19.8407109
61 0.3517835 -10.6702440
62 7.5339502 0.3517835
63 12.3731168 7.5339502
64 12.6780312 12.3731168
65 19.3541978 12.6780312
66 23.2271978 19.3541978
67 25.6086978 23.2271978
68 20.5491978 25.6086978
69 20.1086978 20.5491978
70 24.0400312 20.1086978
71 31.2076978 24.0400312
72 36.8401648 31.2076978
> 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/73dvk1291053163.ps",horizontal=F,onefile=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/83dvk1291053163.ps",horizontal=F,onefile=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/93dvk1291053163.ps",horizontal=F,onefile=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/10dmco1291053163.ps",horizontal=F,onefile=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/11hnbb1291053163.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/12knrh1291053163.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/1396ob1291053163.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/142g5w1291053163.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/155ym21291053163.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/161q2t1291053163.tab")
+ }
>
> try(system("convert tmp/163xc1291053163.ps tmp/163xc1291053163.png",intern=TRUE))
character(0)
> try(system("convert tmp/2zcwx1291053163.ps tmp/2zcwx1291053163.png",intern=TRUE))
character(0)
> try(system("convert tmp/3zcwx1291053163.ps tmp/3zcwx1291053163.png",intern=TRUE))
character(0)
> try(system("convert tmp/4zcwx1291053163.ps tmp/4zcwx1291053163.png",intern=TRUE))
character(0)
> try(system("convert tmp/5zcwx1291053163.ps tmp/5zcwx1291053163.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ame01291053163.ps tmp/6ame01291053163.png",intern=TRUE))
character(0)
> try(system("convert tmp/73dvk1291053163.ps tmp/73dvk1291053163.png",intern=TRUE))
character(0)
> try(system("convert tmp/83dvk1291053163.ps tmp/83dvk1291053163.png",intern=TRUE))
character(0)
> try(system("convert tmp/93dvk1291053163.ps tmp/93dvk1291053163.png",intern=TRUE))
character(0)
> try(system("convert tmp/10dmco1291053163.ps tmp/10dmco1291053163.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.011 2.481 4.516