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(376.974,377.632,378.205,370.861,369.167,371.551,382.842,381.903,384.502,392.058,384.359,388.884,386.586,387.495,385.705,378.67,377.367,376.911,389.827,387.82,387.267,380.575,372.402,376.74,377.795,376.126,370.804,367.98,367.866,366.121,379.421,378.519,372.423,355.072,344.693,342.892,344.178,337.606,327.103,323.953,316.532,306.307,327.225,329.573,313.761,307.836,300.074,304.198,306.122,300.414,292.133,290.616,280.244,285.179,305.486,305.957,293.886,289.441,288.776,299.149,306.532,309.914,313.468,314.901,309.16,316.15,336.544,339.196,326.738,320.838,318.62,331.533,335.378),dim=c(1,73),dimnames=list(c('Maandelijksewerkloosheid'),1:73))
> y <- array(NA,dim=c(1,73),dimnames=list(c('Maandelijksewerkloosheid'),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 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 376.974 1 0 0 0 0 0 0 0 0 0 0 1
2 377.632 0 1 0 0 0 0 0 0 0 0 0 2
3 378.205 0 0 1 0 0 0 0 0 0 0 0 3
4 370.861 0 0 0 1 0 0 0 0 0 0 0 4
5 369.167 0 0 0 0 1 0 0 0 0 0 0 5
6 371.551 0 0 0 0 0 1 0 0 0 0 0 6
7 382.842 0 0 0 0 0 0 1 0 0 0 0 7
8 381.903 0 0 0 0 0 0 0 1 0 0 0 8
9 384.502 0 0 0 0 0 0 0 0 1 0 0 9
10 392.058 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 1 11
12 388.884 0 0 0 0 0 0 0 0 0 0 0 12
13 386.586 1 0 0 0 0 0 0 0 0 0 0 13
14 387.495 0 1 0 0 0 0 0 0 0 0 0 14
15 385.705 0 0 1 0 0 0 0 0 0 0 0 15
16 378.670 0 0 0 1 0 0 0 0 0 0 0 16
17 377.367 0 0 0 0 1 0 0 0 0 0 0 17
18 376.911 0 0 0 0 0 1 0 0 0 0 0 18
19 389.827 0 0 0 0 0 0 1 0 0 0 0 19
20 387.820 0 0 0 0 0 0 0 1 0 0 0 20
21 387.267 0 0 0 0 0 0 0 0 1 0 0 21
22 380.575 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 1 23
24 376.740 0 0 0 0 0 0 0 0 0 0 0 24
25 377.795 1 0 0 0 0 0 0 0 0 0 0 25
26 376.126 0 1 0 0 0 0 0 0 0 0 0 26
27 370.804 0 0 1 0 0 0 0 0 0 0 0 27
28 367.980 0 0 0 1 0 0 0 0 0 0 0 28
29 367.866 0 0 0 0 1 0 0 0 0 0 0 29
30 366.121 0 0 0 0 0 1 0 0 0 0 0 30
31 379.421 0 0 0 0 0 0 1 0 0 0 0 31
32 378.519 0 0 0 0 0 0 0 1 0 0 0 32
33 372.423 0 0 0 0 0 0 0 0 1 0 0 33
34 355.072 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 1 35
36 342.892 0 0 0 0 0 0 0 0 0 0 0 36
37 344.178 1 0 0 0 0 0 0 0 0 0 0 37
38 337.606 0 1 0 0 0 0 0 0 0 0 0 38
39 327.103 0 0 1 0 0 0 0 0 0 0 0 39
40 323.953 0 0 0 1 0 0 0 0 0 0 0 40
41 316.532 0 0 0 0 1 0 0 0 0 0 0 41
42 306.307 0 0 0 0 0 1 0 0 0 0 0 42
43 327.225 0 0 0 0 0 0 1 0 0 0 0 43
44 329.573 0 0 0 0 0 0 0 1 0 0 0 44
45 313.761 0 0 0 0 0 0 0 0 1 0 0 45
46 307.836 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 1 47
48 304.198 0 0 0 0 0 0 0 0 0 0 0 48
49 306.122 1 0 0 0 0 0 0 0 0 0 0 49
50 300.414 0 1 0 0 0 0 0 0 0 0 0 50
51 292.133 0 0 1 0 0 0 0 0 0 0 0 51
52 290.616 0 0 0 1 0 0 0 0 0 0 0 52
53 280.244 0 0 0 0 1 0 0 0 0 0 0 53
54 285.179 0 0 0 0 0 1 0 0 0 0 0 54
55 305.486 0 0 0 0 0 0 1 0 0 0 0 55
56 305.957 0 0 0 0 0 0 0 1 0 0 0 56
57 293.886 0 0 0 0 0 0 0 0 1 0 0 57
58 289.441 0 0 0 0 0 0 0 0 0 1 0 58
59 288.776 0 0 0 0 0 0 0 0 0 0 1 59
60 299.149 0 0 0 0 0 0 0 0 0 0 0 60
61 306.532 1 0 0 0 0 0 0 0 0 0 0 61
62 309.914 0 1 0 0 0 0 0 0 0 0 0 62
63 313.468 0 0 1 0 0 0 0 0 0 0 0 63
64 314.901 0 0 0 1 0 0 0 0 0 0 0 64
65 309.160 0 0 0 0 1 0 0 0 0 0 0 65
66 316.150 0 0 0 0 0 1 0 0 0 0 0 66
67 336.544 0 0 0 0 0 0 1 0 0 0 0 67
68 339.196 0 0 0 0 0 0 0 1 0 0 0 68
69 326.738 0 0 0 0 0 0 0 0 1 0 0 69
70 320.838 0 0 0 0 0 0 0 0 0 1 0 70
71 318.620 0 0 0 0 0 0 0 0 0 0 1 71
72 331.533 0 0 0 0 0 0 0 0 0 0 0 72
73 335.378 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) M1 M2 M3 M4 M5
399.0889 0.1191 -6.3022 -8.5370 -10.5497 -13.5972
M6 M7 M8 M9 M10 M11
-11.8899 6.0245 7.6884 1.6833 -2.3828 -7.1387
t
-1.3934
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-31.397 -16.336 3.518 15.563 37.888
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 399.0889 9.5696 41.704 <2e-16 ***
M1 0.1191 11.3084 0.011 0.992
M2 -6.3022 11.7756 -0.535 0.594
M3 -8.5370 11.7651 -0.726 0.471
M4 -10.5497 11.7557 -0.897 0.373
M5 -13.5972 11.7475 -1.157 0.252
M6 -11.8899 11.7403 -1.013 0.315
M7 6.0245 11.7342 0.513 0.610
M8 7.6884 11.7292 0.655 0.515
M9 1.6833 11.7253 0.144 0.886
M10 -2.3828 11.7226 -0.203 0.840
M11 -7.1387 11.7209 -0.609 0.545
t -1.3934 0.1139 -12.231 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.3 on 60 degrees of freedom
Multiple R-squared: 0.7223, Adjusted R-squared: 0.6667
F-statistic: 13 on 12 and 60 DF, p-value: 1.385e-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,] 7.827163e-05 1.565433e-04 9.999217e-01
[2,] 2.647398e-06 5.294797e-06 9.999974e-01
[3,] 1.072244e-06 2.144487e-06 9.999989e-01
[4,] 5.829486e-08 1.165897e-07 9.999999e-01
[5,] 5.177419e-09 1.035484e-08 1.000000e+00
[6,] 5.375063e-09 1.075013e-08 1.000000e+00
[7,] 5.835560e-06 1.167112e-05 9.999942e-01
[8,] 1.617818e-05 3.235636e-05 9.999838e-01
[9,] 1.883429e-05 3.766858e-05 9.999812e-01
[10,] 7.031284e-06 1.406257e-05 9.999930e-01
[11,] 3.450558e-06 6.901116e-06 9.999965e-01
[12,] 2.823934e-06 5.647869e-06 9.999972e-01
[13,] 1.204257e-06 2.408514e-06 9.999988e-01
[14,] 5.968956e-07 1.193791e-06 9.999994e-01
[15,] 3.911251e-07 7.822503e-07 9.999996e-01
[16,] 2.170721e-07 4.341442e-07 9.999998e-01
[17,] 1.268169e-07 2.536338e-07 9.999999e-01
[18,] 3.085618e-07 6.171236e-07 9.999997e-01
[19,] 2.118616e-05 4.237233e-05 9.999788e-01
[20,] 3.584877e-04 7.169753e-04 9.996415e-01
[21,] 3.569081e-03 7.138163e-03 9.964309e-01
[22,] 9.762727e-03 1.952545e-02 9.902373e-01
[23,] 3.683390e-02 7.366780e-02 9.631661e-01
[24,] 1.158733e-01 2.317465e-01 8.841267e-01
[25,] 2.189340e-01 4.378680e-01 7.810660e-01
[26,] 4.344385e-01 8.688770e-01 5.655615e-01
[27,] 6.043493e-01 7.913014e-01 3.956507e-01
[28,] 6.938891e-01 6.122219e-01 3.061109e-01
[29,] 7.842034e-01 4.315932e-01 2.157966e-01
[30,] 8.919267e-01 2.161465e-01 1.080733e-01
[31,] 9.654619e-01 6.907629e-02 3.453814e-02
[32,] 9.914148e-01 1.717038e-02 8.585192e-03
[33,] 9.981517e-01 3.696525e-03 1.848262e-03
[34,] 9.999499e-01 1.001672e-04 5.008360e-05
[35,] 9.999997e-01 6.508341e-07 3.254170e-07
[36,] 9.999999e-01 1.150174e-07 5.750868e-08
[37,] 1.000000e+00 1.305802e-08 6.529012e-09
[38,] 1.000000e+00 6.485869e-08 3.242935e-08
[39,] 9.999995e-01 1.004586e-06 5.022930e-07
[40,] 9.999925e-01 1.506995e-05 7.534977e-06
[41,] 9.999412e-01 1.175409e-04 5.877044e-05
[42,] 9.995533e-01 8.934082e-04 4.467041e-04
> postscript(file="/var/www/html/rcomp/tmp/1dlht1290975166.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/26dyw1290975166.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/36dyw1290975166.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/46dyw1290975166.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/56dyw1290975166.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 = 73
Frequency = 1
1 2 3 4 5 6 7
-20.840667 -12.367937 -8.166770 -12.104603 -9.357770 -7.287603 -12.517603
8 9 10 11 12 13 14
-13.727103 -3.729603 9.285897 7.736230 6.515897 5.492175 14.215905
15 16 17 18 19 20 21
16.054071 12.425238 15.563071 14.793238 11.188238 8.910738 15.756238
22 23 24 25 26 27 28
14.523738 12.500071 11.092738 13.422016 19.567746 17.873913 18.456079
29 30 31 32 33 34 35
22.782913 20.724079 17.503079 16.330579 17.633079 5.741579 1.511913
36 37 38 39 40 41 42
-6.034421 -3.474143 -2.231413 -9.106246 -8.850079 -11.830246 -22.369079
43 44 45 46 47 48 49
-17.972079 -15.894579 -24.308079 -24.773579 -26.386246 -28.007579 -24.809302
50 51 52 53 54 55 56
-22.702571 -27.355405 -25.466238 -31.397405 -26.776238 -22.990238 -22.789738
57 58 59 60 61 62 63
-27.462238 -26.447738 -20.963405 -16.335738 -7.678460 3.518270 10.700437
64 65 66 67 68 69 70
15.539603 14.239437 20.915603 24.788603 27.170103 22.110603 21.670103
71 72 73
25.601437 32.769103 37.888381
> postscript(file="/var/www/html/rcomp/tmp/6gmfz1290975166.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 = 73
Frequency = 1
lag(myerror, k = 1) myerror
0 -20.840667 NA
1 -12.367937 -20.840667
2 -8.166770 -12.367937
3 -12.104603 -8.166770
4 -9.357770 -12.104603
5 -7.287603 -9.357770
6 -12.517603 -7.287603
7 -13.727103 -12.517603
8 -3.729603 -13.727103
9 9.285897 -3.729603
10 7.736230 9.285897
11 6.515897 7.736230
12 5.492175 6.515897
13 14.215905 5.492175
14 16.054071 14.215905
15 12.425238 16.054071
16 15.563071 12.425238
17 14.793238 15.563071
18 11.188238 14.793238
19 8.910738 11.188238
20 15.756238 8.910738
21 14.523738 15.756238
22 12.500071 14.523738
23 11.092738 12.500071
24 13.422016 11.092738
25 19.567746 13.422016
26 17.873913 19.567746
27 18.456079 17.873913
28 22.782913 18.456079
29 20.724079 22.782913
30 17.503079 20.724079
31 16.330579 17.503079
32 17.633079 16.330579
33 5.741579 17.633079
34 1.511913 5.741579
35 -6.034421 1.511913
36 -3.474143 -6.034421
37 -2.231413 -3.474143
38 -9.106246 -2.231413
39 -8.850079 -9.106246
40 -11.830246 -8.850079
41 -22.369079 -11.830246
42 -17.972079 -22.369079
43 -15.894579 -17.972079
44 -24.308079 -15.894579
45 -24.773579 -24.308079
46 -26.386246 -24.773579
47 -28.007579 -26.386246
48 -24.809302 -28.007579
49 -22.702571 -24.809302
50 -27.355405 -22.702571
51 -25.466238 -27.355405
52 -31.397405 -25.466238
53 -26.776238 -31.397405
54 -22.990238 -26.776238
55 -22.789738 -22.990238
56 -27.462238 -22.789738
57 -26.447738 -27.462238
58 -20.963405 -26.447738
59 -16.335738 -20.963405
60 -7.678460 -16.335738
61 3.518270 -7.678460
62 10.700437 3.518270
63 15.539603 10.700437
64 14.239437 15.539603
65 20.915603 14.239437
66 24.788603 20.915603
67 27.170103 24.788603
68 22.110603 27.170103
69 21.670103 22.110603
70 25.601437 21.670103
71 32.769103 25.601437
72 37.888381 32.769103
73 NA 37.888381
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -12.367937 -20.840667
[2,] -8.166770 -12.367937
[3,] -12.104603 -8.166770
[4,] -9.357770 -12.104603
[5,] -7.287603 -9.357770
[6,] -12.517603 -7.287603
[7,] -13.727103 -12.517603
[8,] -3.729603 -13.727103
[9,] 9.285897 -3.729603
[10,] 7.736230 9.285897
[11,] 6.515897 7.736230
[12,] 5.492175 6.515897
[13,] 14.215905 5.492175
[14,] 16.054071 14.215905
[15,] 12.425238 16.054071
[16,] 15.563071 12.425238
[17,] 14.793238 15.563071
[18,] 11.188238 14.793238
[19,] 8.910738 11.188238
[20,] 15.756238 8.910738
[21,] 14.523738 15.756238
[22,] 12.500071 14.523738
[23,] 11.092738 12.500071
[24,] 13.422016 11.092738
[25,] 19.567746 13.422016
[26,] 17.873913 19.567746
[27,] 18.456079 17.873913
[28,] 22.782913 18.456079
[29,] 20.724079 22.782913
[30,] 17.503079 20.724079
[31,] 16.330579 17.503079
[32,] 17.633079 16.330579
[33,] 5.741579 17.633079
[34,] 1.511913 5.741579
[35,] -6.034421 1.511913
[36,] -3.474143 -6.034421
[37,] -2.231413 -3.474143
[38,] -9.106246 -2.231413
[39,] -8.850079 -9.106246
[40,] -11.830246 -8.850079
[41,] -22.369079 -11.830246
[42,] -17.972079 -22.369079
[43,] -15.894579 -17.972079
[44,] -24.308079 -15.894579
[45,] -24.773579 -24.308079
[46,] -26.386246 -24.773579
[47,] -28.007579 -26.386246
[48,] -24.809302 -28.007579
[49,] -22.702571 -24.809302
[50,] -27.355405 -22.702571
[51,] -25.466238 -27.355405
[52,] -31.397405 -25.466238
[53,] -26.776238 -31.397405
[54,] -22.990238 -26.776238
[55,] -22.789738 -22.990238
[56,] -27.462238 -22.789738
[57,] -26.447738 -27.462238
[58,] -20.963405 -26.447738
[59,] -16.335738 -20.963405
[60,] -7.678460 -16.335738
[61,] 3.518270 -7.678460
[62,] 10.700437 3.518270
[63,] 15.539603 10.700437
[64,] 14.239437 15.539603
[65,] 20.915603 14.239437
[66,] 24.788603 20.915603
[67,] 27.170103 24.788603
[68,] 22.110603 27.170103
[69,] 21.670103 22.110603
[70,] 25.601437 21.670103
[71,] 32.769103 25.601437
[72,] 37.888381 32.769103
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -12.367937 -20.840667
2 -8.166770 -12.367937
3 -12.104603 -8.166770
4 -9.357770 -12.104603
5 -7.287603 -9.357770
6 -12.517603 -7.287603
7 -13.727103 -12.517603
8 -3.729603 -13.727103
9 9.285897 -3.729603
10 7.736230 9.285897
11 6.515897 7.736230
12 5.492175 6.515897
13 14.215905 5.492175
14 16.054071 14.215905
15 12.425238 16.054071
16 15.563071 12.425238
17 14.793238 15.563071
18 11.188238 14.793238
19 8.910738 11.188238
20 15.756238 8.910738
21 14.523738 15.756238
22 12.500071 14.523738
23 11.092738 12.500071
24 13.422016 11.092738
25 19.567746 13.422016
26 17.873913 19.567746
27 18.456079 17.873913
28 22.782913 18.456079
29 20.724079 22.782913
30 17.503079 20.724079
31 16.330579 17.503079
32 17.633079 16.330579
33 5.741579 17.633079
34 1.511913 5.741579
35 -6.034421 1.511913
36 -3.474143 -6.034421
37 -2.231413 -3.474143
38 -9.106246 -2.231413
39 -8.850079 -9.106246
40 -11.830246 -8.850079
41 -22.369079 -11.830246
42 -17.972079 -22.369079
43 -15.894579 -17.972079
44 -24.308079 -15.894579
45 -24.773579 -24.308079
46 -26.386246 -24.773579
47 -28.007579 -26.386246
48 -24.809302 -28.007579
49 -22.702571 -24.809302
50 -27.355405 -22.702571
51 -25.466238 -27.355405
52 -31.397405 -25.466238
53 -26.776238 -31.397405
54 -22.990238 -26.776238
55 -22.789738 -22.990238
56 -27.462238 -22.789738
57 -26.447738 -27.462238
58 -20.963405 -26.447738
59 -16.335738 -20.963405
60 -7.678460 -16.335738
61 3.518270 -7.678460
62 10.700437 3.518270
63 15.539603 10.700437
64 14.239437 15.539603
65 20.915603 14.239437
66 24.788603 20.915603
67 27.170103 24.788603
68 22.110603 27.170103
69 21.670103 22.110603
70 25.601437 21.670103
71 32.769103 25.601437
72 37.888381 32.769103
> 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/7wg4q1290975166.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/8wg4q1290975166.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/9wg4q1290975166.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/102me51290975166.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/1155us1290975166.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/1295tg1290975166.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/13nx8p1290975166.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/148g7v1290975166.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/15uyo11290975166.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/16xhmp1290975166.tab")
+ }
>
> try(system("convert tmp/1dlht1290975166.ps tmp/1dlht1290975166.png",intern=TRUE))
character(0)
> try(system("convert tmp/26dyw1290975166.ps tmp/26dyw1290975166.png",intern=TRUE))
character(0)
> try(system("convert tmp/36dyw1290975166.ps tmp/36dyw1290975166.png",intern=TRUE))
character(0)
> try(system("convert tmp/46dyw1290975166.ps tmp/46dyw1290975166.png",intern=TRUE))
character(0)
> try(system("convert tmp/56dyw1290975166.ps tmp/56dyw1290975166.png",intern=TRUE))
character(0)
> try(system("convert tmp/6gmfz1290975166.ps tmp/6gmfz1290975166.png",intern=TRUE))
character(0)
> try(system("convert tmp/7wg4q1290975166.ps tmp/7wg4q1290975166.png",intern=TRUE))
character(0)
> try(system("convert tmp/8wg4q1290975166.ps tmp/8wg4q1290975166.png",intern=TRUE))
character(0)
> try(system("convert tmp/9wg4q1290975166.ps tmp/9wg4q1290975166.png",intern=TRUE))
character(0)
> try(system("convert tmp/102me51290975166.ps tmp/102me51290975166.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.576 1.574 6.137