R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-redhat-linux-gnu (64-bit)
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(16,17,23,24,27,31,40,47,43,60,64,65,65,55,57,57,57,65,69,70,71,71,73,68,65,57,41,21,21,17,9,11,6,-2,0,5,3,7,4,8,9,14,12,12,7,15,14,19,39,12,11,17,16,25,24,28,25,31,24,24),dim=c(1,60),dimnames=list(c('Werkloosheid'),1:60))
> y <- array(NA,dim=c(1,60),dimnames=list(c('Werkloosheid'),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
> 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
Werkloosheid M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 16 1 0 0 0 0 0 0 0 0 0 0 1
2 17 0 1 0 0 0 0 0 0 0 0 0 2
3 23 0 0 1 0 0 0 0 0 0 0 0 3
4 24 0 0 0 1 0 0 0 0 0 0 0 4
5 27 0 0 0 0 1 0 0 0 0 0 0 5
6 31 0 0 0 0 0 1 0 0 0 0 0 6
7 40 0 0 0 0 0 0 1 0 0 0 0 7
8 47 0 0 0 0 0 0 0 1 0 0 0 8
9 43 0 0 0 0 0 0 0 0 1 0 0 9
10 60 0 0 0 0 0 0 0 0 0 1 0 10
11 64 0 0 0 0 0 0 0 0 0 0 1 11
12 65 0 0 0 0 0 0 0 0 0 0 0 12
13 65 1 0 0 0 0 0 0 0 0 0 0 13
14 55 0 1 0 0 0 0 0 0 0 0 0 14
15 57 0 0 1 0 0 0 0 0 0 0 0 15
16 57 0 0 0 1 0 0 0 0 0 0 0 16
17 57 0 0 0 0 1 0 0 0 0 0 0 17
18 65 0 0 0 0 0 1 0 0 0 0 0 18
19 69 0 0 0 0 0 0 1 0 0 0 0 19
20 70 0 0 0 0 0 0 0 1 0 0 0 20
21 71 0 0 0 0 0 0 0 0 1 0 0 21
22 71 0 0 0 0 0 0 0 0 0 1 0 22
23 73 0 0 0 0 0 0 0 0 0 0 1 23
24 68 0 0 0 0 0 0 0 0 0 0 0 24
25 65 1 0 0 0 0 0 0 0 0 0 0 25
26 57 0 1 0 0 0 0 0 0 0 0 0 26
27 41 0 0 1 0 0 0 0 0 0 0 0 27
28 21 0 0 0 1 0 0 0 0 0 0 0 28
29 21 0 0 0 0 1 0 0 0 0 0 0 29
30 17 0 0 0 0 0 1 0 0 0 0 0 30
31 9 0 0 0 0 0 0 1 0 0 0 0 31
32 11 0 0 0 0 0 0 0 1 0 0 0 32
33 6 0 0 0 0 0 0 0 0 1 0 0 33
34 -2 0 0 0 0 0 0 0 0 0 1 0 34
35 0 0 0 0 0 0 0 0 0 0 0 1 35
36 5 0 0 0 0 0 0 0 0 0 0 0 36
37 3 1 0 0 0 0 0 0 0 0 0 0 37
38 7 0 1 0 0 0 0 0 0 0 0 0 38
39 4 0 0 1 0 0 0 0 0 0 0 0 39
40 8 0 0 0 1 0 0 0 0 0 0 0 40
41 9 0 0 0 0 1 0 0 0 0 0 0 41
42 14 0 0 0 0 0 1 0 0 0 0 0 42
43 12 0 0 0 0 0 0 1 0 0 0 0 43
44 12 0 0 0 0 0 0 0 1 0 0 0 44
45 7 0 0 0 0 0 0 0 0 1 0 0 45
46 15 0 0 0 0 0 0 0 0 0 1 0 46
47 14 0 0 0 0 0 0 0 0 0 0 1 47
48 19 0 0 0 0 0 0 0 0 0 0 0 48
49 39 1 0 0 0 0 0 0 0 0 0 0 49
50 12 0 1 0 0 0 0 0 0 0 0 0 50
51 11 0 0 1 0 0 0 0 0 0 0 0 51
52 17 0 0 0 1 0 0 0 0 0 0 0 52
53 16 0 0 0 0 1 0 0 0 0 0 0 53
54 25 0 0 0 0 0 1 0 0 0 0 0 54
55 24 0 0 0 0 0 0 1 0 0 0 0 55
56 28 0 0 0 0 0 0 0 1 0 0 0 56
57 25 0 0 0 0 0 0 0 0 1 0 0 57
58 31 0 0 0 0 0 0 0 0 0 1 0 58
59 24 0 0 0 0 0 0 0 0 0 0 1 59
60 24 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) M1 M2 M3 M4 M5
61.6000 -6.3611 -13.6556 -15.3500 -16.4444 -15.1389
M6 M7 M8 M9 M10 M11
-10.0333 -8.9278 -5.4222 -7.9167 -2.6111 -1.9056
t
-0.7056
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-38.53 -14.28 -2.10 14.58 32.13
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 61.6000 11.4046 5.401 2.14e-06 ***
M1 -6.3611 13.8744 -0.458 0.649
M2 -13.6556 13.8537 -0.986 0.329
M3 -15.3500 13.8349 -1.110 0.273
M4 -16.4444 13.8181 -1.190 0.240
M5 -15.1389 13.8032 -1.097 0.278
M6 -10.0333 13.7903 -0.728 0.470
M7 -8.9278 13.7794 -0.648 0.520
M8 -5.4222 13.7704 -0.394 0.696
M9 -7.9167 13.7635 -0.575 0.568
M10 -2.6111 13.7585 -0.190 0.850
M11 -1.9056 13.7555 -0.139 0.890
t -0.7056 0.1654 -4.265 9.59e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.75 on 47 degrees of freedom
Multiple R-squared: 0.2995, Adjusted R-squared: 0.1207
F-statistic: 1.675 on 12 and 47 DF, p-value: 0.1036
> 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.0171461857 3.429237e-02 9.828538e-01
[2,] 0.0058794922 1.175898e-02 9.941205e-01
[3,] 0.0012975051 2.595010e-03 9.987025e-01
[4,] 0.0004687397 9.374794e-04 9.995313e-01
[5,] 0.0004053696 8.107392e-04 9.995946e-01
[6,] 0.0001726824 3.453649e-04 9.998273e-01
[7,] 0.0011842805 2.368561e-03 9.988157e-01
[8,] 0.0054676046 1.093521e-02 9.945324e-01
[9,] 0.0318365759 6.367315e-02 9.681634e-01
[10,] 0.0732417805 1.464836e-01 9.267582e-01
[11,] 0.3219506784 6.439014e-01 6.780493e-01
[12,] 0.9262217843 1.475564e-01 7.377822e-02
[13,] 0.9973761987 5.247603e-03 2.623801e-03
[14,] 0.9999330938 1.338125e-04 6.690623e-05
[15,] 0.9999926139 1.477229e-05 7.386146e-06
[16,] 0.9999981944 3.611238e-06 1.805619e-06
[17,] 0.9999991598 1.680382e-06 8.401910e-07
[18,] 0.9999994344 1.131256e-06 5.656281e-07
[19,] 0.9999994913 1.017408e-06 5.087042e-07
[20,] 0.9999989777 2.044571e-06 1.022286e-06
[21,] 0.9999970149 5.970257e-06 2.985128e-06
[22,] 0.9999999818 3.645068e-08 1.822534e-08
[23,] 0.9999999274 1.451288e-07 7.256441e-08
[24,] 0.9999995803 8.393420e-07 4.196710e-07
[25,] 0.9999967376 6.524765e-06 3.262383e-06
[26,] 0.9999845055 3.098902e-05 1.549451e-05
[27,] 0.9998714231 2.571538e-04 1.285769e-04
[28,] 0.9989780469 2.043906e-03 1.021953e-03
[29,] 0.9936698515 1.266030e-02 6.330149e-03
> postscript(file="/var/www/wessaorg/rcomp/tmp/1d1g11293648893.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/wessaorg/rcomp/tmp/2obxm1293648893.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/wessaorg/rcomp/tmp/3obxm1293648893.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/wessaorg/rcomp/tmp/4obxm1293648893.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/wessaorg/rcomp/tmp/5obxm1293648893.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 = 60
Frequency = 1
1 2 3 4 5 6
-38.5333333 -29.5333333 -21.1333333 -18.3333333 -15.9333333 -16.3333333
7 8 9 10 11 12
-7.7333333 -3.5333333 -4.3333333 8.0666667 12.0666667 11.8666667
13 14 15 16 17 18
18.9333333 16.9333333 21.3333333 23.1333333 22.5333333 26.1333333
19 20 21 22 23 24
29.7333333 27.9333333 32.1333333 27.5333333 29.5333333 23.3333333
25 26 27 28 29 30
27.4000000 27.4000000 13.8000000 -4.4000000 -5.0000000 -13.4000000
31 32 33 34 35 36
-21.8000000 -22.6000000 -24.4000000 -37.0000000 -35.0000000 -31.2000000
37 38 39 40 41 42
-26.1333333 -14.1333333 -14.7333333 -8.9333333 -8.5333333 -7.9333333
43 44 45 46 47 48
-10.3333333 -13.1333333 -14.9333333 -11.5333333 -12.5333333 -8.7333333
49 50 51 52 53 54
18.3333333 -0.6666667 0.7333333 8.5333333 6.9333333 11.5333333
55 56 57 58 59 60
10.1333333 11.3333333 11.5333333 12.9333333 5.9333333 4.7333333
> postscript(file="/var/www/wessaorg/rcomp/tmp/6gke71293648893.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 -38.5333333 NA
1 -29.5333333 -38.5333333
2 -21.1333333 -29.5333333
3 -18.3333333 -21.1333333
4 -15.9333333 -18.3333333
5 -16.3333333 -15.9333333
6 -7.7333333 -16.3333333
7 -3.5333333 -7.7333333
8 -4.3333333 -3.5333333
9 8.0666667 -4.3333333
10 12.0666667 8.0666667
11 11.8666667 12.0666667
12 18.9333333 11.8666667
13 16.9333333 18.9333333
14 21.3333333 16.9333333
15 23.1333333 21.3333333
16 22.5333333 23.1333333
17 26.1333333 22.5333333
18 29.7333333 26.1333333
19 27.9333333 29.7333333
20 32.1333333 27.9333333
21 27.5333333 32.1333333
22 29.5333333 27.5333333
23 23.3333333 29.5333333
24 27.4000000 23.3333333
25 27.4000000 27.4000000
26 13.8000000 27.4000000
27 -4.4000000 13.8000000
28 -5.0000000 -4.4000000
29 -13.4000000 -5.0000000
30 -21.8000000 -13.4000000
31 -22.6000000 -21.8000000
32 -24.4000000 -22.6000000
33 -37.0000000 -24.4000000
34 -35.0000000 -37.0000000
35 -31.2000000 -35.0000000
36 -26.1333333 -31.2000000
37 -14.1333333 -26.1333333
38 -14.7333333 -14.1333333
39 -8.9333333 -14.7333333
40 -8.5333333 -8.9333333
41 -7.9333333 -8.5333333
42 -10.3333333 -7.9333333
43 -13.1333333 -10.3333333
44 -14.9333333 -13.1333333
45 -11.5333333 -14.9333333
46 -12.5333333 -11.5333333
47 -8.7333333 -12.5333333
48 18.3333333 -8.7333333
49 -0.6666667 18.3333333
50 0.7333333 -0.6666667
51 8.5333333 0.7333333
52 6.9333333 8.5333333
53 11.5333333 6.9333333
54 10.1333333 11.5333333
55 11.3333333 10.1333333
56 11.5333333 11.3333333
57 12.9333333 11.5333333
58 5.9333333 12.9333333
59 4.7333333 5.9333333
60 NA 4.7333333
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -29.5333333 -38.5333333
[2,] -21.1333333 -29.5333333
[3,] -18.3333333 -21.1333333
[4,] -15.9333333 -18.3333333
[5,] -16.3333333 -15.9333333
[6,] -7.7333333 -16.3333333
[7,] -3.5333333 -7.7333333
[8,] -4.3333333 -3.5333333
[9,] 8.0666667 -4.3333333
[10,] 12.0666667 8.0666667
[11,] 11.8666667 12.0666667
[12,] 18.9333333 11.8666667
[13,] 16.9333333 18.9333333
[14,] 21.3333333 16.9333333
[15,] 23.1333333 21.3333333
[16,] 22.5333333 23.1333333
[17,] 26.1333333 22.5333333
[18,] 29.7333333 26.1333333
[19,] 27.9333333 29.7333333
[20,] 32.1333333 27.9333333
[21,] 27.5333333 32.1333333
[22,] 29.5333333 27.5333333
[23,] 23.3333333 29.5333333
[24,] 27.4000000 23.3333333
[25,] 27.4000000 27.4000000
[26,] 13.8000000 27.4000000
[27,] -4.4000000 13.8000000
[28,] -5.0000000 -4.4000000
[29,] -13.4000000 -5.0000000
[30,] -21.8000000 -13.4000000
[31,] -22.6000000 -21.8000000
[32,] -24.4000000 -22.6000000
[33,] -37.0000000 -24.4000000
[34,] -35.0000000 -37.0000000
[35,] -31.2000000 -35.0000000
[36,] -26.1333333 -31.2000000
[37,] -14.1333333 -26.1333333
[38,] -14.7333333 -14.1333333
[39,] -8.9333333 -14.7333333
[40,] -8.5333333 -8.9333333
[41,] -7.9333333 -8.5333333
[42,] -10.3333333 -7.9333333
[43,] -13.1333333 -10.3333333
[44,] -14.9333333 -13.1333333
[45,] -11.5333333 -14.9333333
[46,] -12.5333333 -11.5333333
[47,] -8.7333333 -12.5333333
[48,] 18.3333333 -8.7333333
[49,] -0.6666667 18.3333333
[50,] 0.7333333 -0.6666667
[51,] 8.5333333 0.7333333
[52,] 6.9333333 8.5333333
[53,] 11.5333333 6.9333333
[54,] 10.1333333 11.5333333
[55,] 11.3333333 10.1333333
[56,] 11.5333333 11.3333333
[57,] 12.9333333 11.5333333
[58,] 5.9333333 12.9333333
[59,] 4.7333333 5.9333333
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -29.5333333 -38.5333333
2 -21.1333333 -29.5333333
3 -18.3333333 -21.1333333
4 -15.9333333 -18.3333333
5 -16.3333333 -15.9333333
6 -7.7333333 -16.3333333
7 -3.5333333 -7.7333333
8 -4.3333333 -3.5333333
9 8.0666667 -4.3333333
10 12.0666667 8.0666667
11 11.8666667 12.0666667
12 18.9333333 11.8666667
13 16.9333333 18.9333333
14 21.3333333 16.9333333
15 23.1333333 21.3333333
16 22.5333333 23.1333333
17 26.1333333 22.5333333
18 29.7333333 26.1333333
19 27.9333333 29.7333333
20 32.1333333 27.9333333
21 27.5333333 32.1333333
22 29.5333333 27.5333333
23 23.3333333 29.5333333
24 27.4000000 23.3333333
25 27.4000000 27.4000000
26 13.8000000 27.4000000
27 -4.4000000 13.8000000
28 -5.0000000 -4.4000000
29 -13.4000000 -5.0000000
30 -21.8000000 -13.4000000
31 -22.6000000 -21.8000000
32 -24.4000000 -22.6000000
33 -37.0000000 -24.4000000
34 -35.0000000 -37.0000000
35 -31.2000000 -35.0000000
36 -26.1333333 -31.2000000
37 -14.1333333 -26.1333333
38 -14.7333333 -14.1333333
39 -8.9333333 -14.7333333
40 -8.5333333 -8.9333333
41 -7.9333333 -8.5333333
42 -10.3333333 -7.9333333
43 -13.1333333 -10.3333333
44 -14.9333333 -13.1333333
45 -11.5333333 -14.9333333
46 -12.5333333 -11.5333333
47 -8.7333333 -12.5333333
48 18.3333333 -8.7333333
49 -0.6666667 18.3333333
50 0.7333333 -0.6666667
51 8.5333333 0.7333333
52 6.9333333 8.5333333
53 11.5333333 6.9333333
54 10.1333333 11.5333333
55 11.3333333 10.1333333
56 11.5333333 11.3333333
57 12.9333333 11.5333333
58 5.9333333 12.9333333
59 4.7333333 5.9333333
> 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/wessaorg/rcomp/tmp/7rtes1293648893.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/wessaorg/rcomp/tmp/8rtes1293648893.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/wessaorg/rcomp/tmp/9rtes1293648893.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/wessaorg/rcomp/tmp/1023dv1293648893.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/wessaorg/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/wessaorg/rcomp/tmp/115lt11293648893.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/wessaorg/rcomp/tmp/12q3s71293648893.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/wessaorg/rcomp/tmp/135vqg1293648893.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/wessaorg/rcomp/tmp/148wom1293648893.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/wessaorg/rcomp/tmp/15cxna1293648893.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/wessaorg/rcomp/tmp/16fx3f1293648893.tab")
+ }
>
> try(system("convert tmp/1d1g11293648893.ps tmp/1d1g11293648893.png",intern=TRUE))
character(0)
> try(system("convert tmp/2obxm1293648893.ps tmp/2obxm1293648893.png",intern=TRUE))
character(0)
> try(system("convert tmp/3obxm1293648893.ps tmp/3obxm1293648893.png",intern=TRUE))
character(0)
> try(system("convert tmp/4obxm1293648893.ps tmp/4obxm1293648893.png",intern=TRUE))
character(0)
> try(system("convert tmp/5obxm1293648893.ps tmp/5obxm1293648893.png",intern=TRUE))
character(0)
> try(system("convert tmp/6gke71293648893.ps tmp/6gke71293648893.png",intern=TRUE))
character(0)
> try(system("convert tmp/7rtes1293648893.ps tmp/7rtes1293648893.png",intern=TRUE))
character(0)
> try(system("convert tmp/8rtes1293648893.ps tmp/8rtes1293648893.png",intern=TRUE))
character(0)
> try(system("convert tmp/9rtes1293648893.ps tmp/9rtes1293648893.png",intern=TRUE))
character(0)
> try(system("convert tmp/1023dv1293648893.ps tmp/1023dv1293648893.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
2.93 0.50 3.59