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.
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(269645,0,267037,0,258113,0,262813,0,267413,0,267366,0,264777,0,258863,0,254844,0,254868,0,277267,0,285351,0,286602,0,283042,0,276687,0,277915,0,277128,0,277103,0,275037,0,270150,0,267140,0,264993,0,287259,0,291186,0,292300,0,288186,0,281477,0,282656,1,280190,1,280408,1,276836,1,275216,1,274352,1,271311,1,289802,1,290726,1,292300,1,278506,1,269826,1,265861,1,269034,1,264176,1,255198,1,253353,1,246057,1,235372,1,258556,1,260993,1,254663,1,250643,1,243422,1,247105,1,248541,1,245039,1,237080,1,237085,1,225554,1,226839,1,247934,1,248333,1,246969,1,245098,1,246263,1),dim=c(2,63),dimnames=list(c('Mannen','Dummy'),1:63))
> y <- array(NA,dim=c(2,63),dimnames=list(c('Mannen','Dummy'),1:63))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Mannen Dummy M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 269645 0 1 0 0 0 0 0 0 0 0 0 0
2 267037 0 0 1 0 0 0 0 0 0 0 0 0
3 258113 0 0 0 1 0 0 0 0 0 0 0 0
4 262813 0 0 0 0 1 0 0 0 0 0 0 0
5 267413 0 0 0 0 0 1 0 0 0 0 0 0
6 267366 0 0 0 0 0 0 1 0 0 0 0 0
7 264777 0 0 0 0 0 0 0 1 0 0 0 0
8 258863 0 0 0 0 0 0 0 0 1 0 0 0
9 254844 0 0 0 0 0 0 0 0 0 1 0 0
10 254868 0 0 0 0 0 0 0 0 0 0 1 0
11 277267 0 0 0 0 0 0 0 0 0 0 0 1
12 285351 0 0 0 0 0 0 0 0 0 0 0 0
13 286602 0 1 0 0 0 0 0 0 0 0 0 0
14 283042 0 0 1 0 0 0 0 0 0 0 0 0
15 276687 0 0 0 1 0 0 0 0 0 0 0 0
16 277915 0 0 0 0 1 0 0 0 0 0 0 0
17 277128 0 0 0 0 0 1 0 0 0 0 0 0
18 277103 0 0 0 0 0 0 1 0 0 0 0 0
19 275037 0 0 0 0 0 0 0 1 0 0 0 0
20 270150 0 0 0 0 0 0 0 0 1 0 0 0
21 267140 0 0 0 0 0 0 0 0 0 1 0 0
22 264993 0 0 0 0 0 0 0 0 0 0 1 0
23 287259 0 0 0 0 0 0 0 0 0 0 0 1
24 291186 0 0 0 0 0 0 0 0 0 0 0 0
25 292300 0 1 0 0 0 0 0 0 0 0 0 0
26 288186 0 0 1 0 0 0 0 0 0 0 0 0
27 281477 0 0 0 1 0 0 0 0 0 0 0 0
28 282656 1 0 0 0 1 0 0 0 0 0 0 0
29 280190 1 0 0 0 0 1 0 0 0 0 0 0
30 280408 1 0 0 0 0 0 1 0 0 0 0 0
31 276836 1 0 0 0 0 0 0 1 0 0 0 0
32 275216 1 0 0 0 0 0 0 0 1 0 0 0
33 274352 1 0 0 0 0 0 0 0 0 1 0 0
34 271311 1 0 0 0 0 0 0 0 0 0 1 0
35 289802 1 0 0 0 0 0 0 0 0 0 0 1
36 290726 1 0 0 0 0 0 0 0 0 0 0 0
37 292300 1 1 0 0 0 0 0 0 0 0 0 0
38 278506 1 0 1 0 0 0 0 0 0 0 0 0
39 269826 1 0 0 1 0 0 0 0 0 0 0 0
40 265861 1 0 0 0 1 0 0 0 0 0 0 0
41 269034 1 0 0 0 0 1 0 0 0 0 0 0
42 264176 1 0 0 0 0 0 1 0 0 0 0 0
43 255198 1 0 0 0 0 0 0 1 0 0 0 0
44 253353 1 0 0 0 0 0 0 0 1 0 0 0
45 246057 1 0 0 0 0 0 0 0 0 1 0 0
46 235372 1 0 0 0 0 0 0 0 0 0 1 0
47 258556 1 0 0 0 0 0 0 0 0 0 0 1
48 260993 1 0 0 0 0 0 0 0 0 0 0 0
49 254663 1 1 0 0 0 0 0 0 0 0 0 0
50 250643 1 0 1 0 0 0 0 0 0 0 0 0
51 243422 1 0 0 1 0 0 0 0 0 0 0 0
52 247105 1 0 0 0 1 0 0 0 0 0 0 0
53 248541 1 0 0 0 0 1 0 0 0 0 0 0
54 245039 1 0 0 0 0 0 1 0 0 0 0 0
55 237080 1 0 0 0 0 0 0 1 0 0 0 0
56 237085 1 0 0 0 0 0 0 0 1 0 0 0
57 225554 1 0 0 0 0 0 0 0 0 1 0 0
58 226839 1 0 0 0 0 0 0 0 0 0 1 0
59 247934 1 0 0 0 0 0 0 0 0 0 0 1
60 248333 1 0 0 0 0 0 0 0 0 0 0 0
61 246969 1 1 0 0 0 0 0 0 0 0 0 0
62 245098 1 0 1 0 0 0 0 0 0 0 0 0
63 246263 1 0 0 1 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Dummy M1 M2 M3 M4
283912 -14323 -3004 -7998 -14119 -8048
M5 M6 M7 M8 M9 M10
-6857 -8499 -13532 -16384 -21728 -24641
M11
-3154
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-22306.24 -11105.25 73.05 7201.50 26491.76
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 283912 7182 39.530 < 2e-16 ***
Dummy -14323 3884 -3.688 0.000557 ***
M1 -3004 9207 -0.326 0.745608
M2 -7998 9207 -0.869 0.389154
M3 -14119 9207 -1.534 0.131453
M4 -8048 9608 -0.838 0.406216
M5 -6857 9608 -0.714 0.478755
M6 -8499 9608 -0.885 0.380580
M7 -13532 9608 -1.408 0.165176
M8 -16384 9608 -1.705 0.094335 .
M9 -21728 9608 -2.262 0.028105 *
M10 -24641 9608 -2.565 0.013374 *
M11 -3154 9608 -0.328 0.744055
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15190 on 50 degrees of freedom
Multiple R-squared: 0.361, Adjusted R-squared: 0.2076
F-statistic: 2.354 on 12 and 50 DF, p-value: 0.01739
> 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.3995664135 0.7991328271 0.6004336
[2,] 0.2652611655 0.5305223310 0.7347388
[3,] 0.1699815345 0.3399630689 0.8300185
[4,] 0.1067271853 0.2134543705 0.8932728
[5,] 0.0686399471 0.1372798942 0.9313601
[6,] 0.0450005793 0.0900011586 0.9549994
[7,] 0.0266355199 0.0532710397 0.9733645
[8,] 0.0153359950 0.0306719899 0.9846640
[9,] 0.0074918013 0.0149836025 0.9925082
[10,] 0.0058128706 0.0116257413 0.9941871
[11,] 0.0040833679 0.0081667358 0.9959166
[12,] 0.0029954511 0.0059909022 0.9970045
[13,] 0.0016720908 0.0033441815 0.9983279
[14,] 0.0008716745 0.0017433490 0.9991283
[15,] 0.0004864482 0.0009728963 0.9995136
[16,] 0.0003220594 0.0006441188 0.9996779
[17,] 0.0002178038 0.0004356077 0.9997822
[18,] 0.0002495189 0.0004990377 0.9997505
[19,] 0.0003811572 0.0007623144 0.9996188
[20,] 0.0006055137 0.0012110274 0.9993945
[21,] 0.0014318160 0.0028636320 0.9985682
[22,] 0.0083232529 0.0166465059 0.9916767
[23,] 0.0296999429 0.0593998858 0.9703001
[24,] 0.0732970198 0.1465940397 0.9267030
[25,] 0.1065861094 0.2131722187 0.8934139
[26,] 0.1560667698 0.3121335396 0.8439332
[27,] 0.2420778518 0.4841557036 0.7579221
[28,] 0.3881392131 0.7762784262 0.6118608
[29,] 0.4912920405 0.9825840811 0.5087080
[30,] 0.7615430568 0.4769138863 0.2384569
[31,] 0.7703250011 0.4593499978 0.2296750
[32,] 0.7715685619 0.4568628763 0.2284314
> postscript(file="/var/www/html/freestat/rcomp/tmp/1hyb71229461042.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/26hni1229461042.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/3k2lf1229461042.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/4j2681229461042.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/5ndmg1229461042.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 = 63
Frequency = 1
1 2 3 4 5 6
-11262.9542 -8876.4542 -11679.7876 -13050.7451 -9641.9451 -8046.1451
7 8 9 10 11 12
-5602.3451 -8664.1451 -7339.1451 -4402.3451 -3490.3451 1439.4549
13 14 15 16 17 18
5694.0458 7128.5458 6894.2124 2051.2549 73.0549 1690.8549
19 20 21 22 23 24
4657.6549 2622.8549 4956.8549 5722.6549 6501.6549 7274.4549
25 26 27 28 29 30
11392.0458 12272.5458 11684.2124 21115.1634 17457.9634 19318.7634
31 32 33 34 35 36
20779.5634 22011.7634 26491.7634 26363.5634 23367.5634 21137.3634
37 38 39 40 41 42
25714.9542 16915.4542 14356.1209 4320.1634 6301.9634 3086.7634
43 44 45 46 47 48
-858.4366 148.7634 -1803.2366 -9575.4366 -7878.4366 -8595.6366
49 50 51 52 53 54
-11922.0458 -10947.5458 -12047.8791 -14435.8366 -14191.0366 -16050.2366
55 56 57 58 59 60
-18976.4366 -16119.2366 -22306.2366 -18108.4366 -18500.4366 -21255.6366
61 62 63
-19616.0458 -16492.5458 -9206.8791
> postscript(file="/var/www/html/freestat/rcomp/tmp/681qv1229461042.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 = 63
Frequency = 1
lag(myerror, k = 1) myerror
0 -11262.9542 NA
1 -8876.4542 -11262.9542
2 -11679.7876 -8876.4542
3 -13050.7451 -11679.7876
4 -9641.9451 -13050.7451
5 -8046.1451 -9641.9451
6 -5602.3451 -8046.1451
7 -8664.1451 -5602.3451
8 -7339.1451 -8664.1451
9 -4402.3451 -7339.1451
10 -3490.3451 -4402.3451
11 1439.4549 -3490.3451
12 5694.0458 1439.4549
13 7128.5458 5694.0458
14 6894.2124 7128.5458
15 2051.2549 6894.2124
16 73.0549 2051.2549
17 1690.8549 73.0549
18 4657.6549 1690.8549
19 2622.8549 4657.6549
20 4956.8549 2622.8549
21 5722.6549 4956.8549
22 6501.6549 5722.6549
23 7274.4549 6501.6549
24 11392.0458 7274.4549
25 12272.5458 11392.0458
26 11684.2124 12272.5458
27 21115.1634 11684.2124
28 17457.9634 21115.1634
29 19318.7634 17457.9634
30 20779.5634 19318.7634
31 22011.7634 20779.5634
32 26491.7634 22011.7634
33 26363.5634 26491.7634
34 23367.5634 26363.5634
35 21137.3634 23367.5634
36 25714.9542 21137.3634
37 16915.4542 25714.9542
38 14356.1209 16915.4542
39 4320.1634 14356.1209
40 6301.9634 4320.1634
41 3086.7634 6301.9634
42 -858.4366 3086.7634
43 148.7634 -858.4366
44 -1803.2366 148.7634
45 -9575.4366 -1803.2366
46 -7878.4366 -9575.4366
47 -8595.6366 -7878.4366
48 -11922.0458 -8595.6366
49 -10947.5458 -11922.0458
50 -12047.8791 -10947.5458
51 -14435.8366 -12047.8791
52 -14191.0366 -14435.8366
53 -16050.2366 -14191.0366
54 -18976.4366 -16050.2366
55 -16119.2366 -18976.4366
56 -22306.2366 -16119.2366
57 -18108.4366 -22306.2366
58 -18500.4366 -18108.4366
59 -21255.6366 -18500.4366
60 -19616.0458 -21255.6366
61 -16492.5458 -19616.0458
62 -9206.8791 -16492.5458
63 NA -9206.8791
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -8876.4542 -11262.9542
[2,] -11679.7876 -8876.4542
[3,] -13050.7451 -11679.7876
[4,] -9641.9451 -13050.7451
[5,] -8046.1451 -9641.9451
[6,] -5602.3451 -8046.1451
[7,] -8664.1451 -5602.3451
[8,] -7339.1451 -8664.1451
[9,] -4402.3451 -7339.1451
[10,] -3490.3451 -4402.3451
[11,] 1439.4549 -3490.3451
[12,] 5694.0458 1439.4549
[13,] 7128.5458 5694.0458
[14,] 6894.2124 7128.5458
[15,] 2051.2549 6894.2124
[16,] 73.0549 2051.2549
[17,] 1690.8549 73.0549
[18,] 4657.6549 1690.8549
[19,] 2622.8549 4657.6549
[20,] 4956.8549 2622.8549
[21,] 5722.6549 4956.8549
[22,] 6501.6549 5722.6549
[23,] 7274.4549 6501.6549
[24,] 11392.0458 7274.4549
[25,] 12272.5458 11392.0458
[26,] 11684.2124 12272.5458
[27,] 21115.1634 11684.2124
[28,] 17457.9634 21115.1634
[29,] 19318.7634 17457.9634
[30,] 20779.5634 19318.7634
[31,] 22011.7634 20779.5634
[32,] 26491.7634 22011.7634
[33,] 26363.5634 26491.7634
[34,] 23367.5634 26363.5634
[35,] 21137.3634 23367.5634
[36,] 25714.9542 21137.3634
[37,] 16915.4542 25714.9542
[38,] 14356.1209 16915.4542
[39,] 4320.1634 14356.1209
[40,] 6301.9634 4320.1634
[41,] 3086.7634 6301.9634
[42,] -858.4366 3086.7634
[43,] 148.7634 -858.4366
[44,] -1803.2366 148.7634
[45,] -9575.4366 -1803.2366
[46,] -7878.4366 -9575.4366
[47,] -8595.6366 -7878.4366
[48,] -11922.0458 -8595.6366
[49,] -10947.5458 -11922.0458
[50,] -12047.8791 -10947.5458
[51,] -14435.8366 -12047.8791
[52,] -14191.0366 -14435.8366
[53,] -16050.2366 -14191.0366
[54,] -18976.4366 -16050.2366
[55,] -16119.2366 -18976.4366
[56,] -22306.2366 -16119.2366
[57,] -18108.4366 -22306.2366
[58,] -18500.4366 -18108.4366
[59,] -21255.6366 -18500.4366
[60,] -19616.0458 -21255.6366
[61,] -16492.5458 -19616.0458
[62,] -9206.8791 -16492.5458
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -8876.4542 -11262.9542
2 -11679.7876 -8876.4542
3 -13050.7451 -11679.7876
4 -9641.9451 -13050.7451
5 -8046.1451 -9641.9451
6 -5602.3451 -8046.1451
7 -8664.1451 -5602.3451
8 -7339.1451 -8664.1451
9 -4402.3451 -7339.1451
10 -3490.3451 -4402.3451
11 1439.4549 -3490.3451
12 5694.0458 1439.4549
13 7128.5458 5694.0458
14 6894.2124 7128.5458
15 2051.2549 6894.2124
16 73.0549 2051.2549
17 1690.8549 73.0549
18 4657.6549 1690.8549
19 2622.8549 4657.6549
20 4956.8549 2622.8549
21 5722.6549 4956.8549
22 6501.6549 5722.6549
23 7274.4549 6501.6549
24 11392.0458 7274.4549
25 12272.5458 11392.0458
26 11684.2124 12272.5458
27 21115.1634 11684.2124
28 17457.9634 21115.1634
29 19318.7634 17457.9634
30 20779.5634 19318.7634
31 22011.7634 20779.5634
32 26491.7634 22011.7634
33 26363.5634 26491.7634
34 23367.5634 26363.5634
35 21137.3634 23367.5634
36 25714.9542 21137.3634
37 16915.4542 25714.9542
38 14356.1209 16915.4542
39 4320.1634 14356.1209
40 6301.9634 4320.1634
41 3086.7634 6301.9634
42 -858.4366 3086.7634
43 148.7634 -858.4366
44 -1803.2366 148.7634
45 -9575.4366 -1803.2366
46 -7878.4366 -9575.4366
47 -8595.6366 -7878.4366
48 -11922.0458 -8595.6366
49 -10947.5458 -11922.0458
50 -12047.8791 -10947.5458
51 -14435.8366 -12047.8791
52 -14191.0366 -14435.8366
53 -16050.2366 -14191.0366
54 -18976.4366 -16050.2366
55 -16119.2366 -18976.4366
56 -22306.2366 -16119.2366
57 -18108.4366 -22306.2366
58 -18500.4366 -18108.4366
59 -21255.6366 -18500.4366
60 -19616.0458 -21255.6366
61 -16492.5458 -19616.0458
62 -9206.8791 -16492.5458
> 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/7hhkl1229461042.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/89o5k1229461042.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/92wq31229461042.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/10wg5q1229461043.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/11gzpv1229461043.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/12az8y1229461043.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/13j8q21229461043.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/14c0zi1229461043.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/15qz8u1229461043.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/16xq2l1229461043.tab")
+ }
>
> system("convert tmp/1hyb71229461042.ps tmp/1hyb71229461042.png")
> system("convert tmp/26hni1229461042.ps tmp/26hni1229461042.png")
> system("convert tmp/3k2lf1229461042.ps tmp/3k2lf1229461042.png")
> system("convert tmp/4j2681229461042.ps tmp/4j2681229461042.png")
> system("convert tmp/5ndmg1229461042.ps tmp/5ndmg1229461042.png")
> system("convert tmp/681qv1229461042.ps tmp/681qv1229461042.png")
> system("convert tmp/7hhkl1229461042.ps tmp/7hhkl1229461042.png")
> system("convert tmp/89o5k1229461042.ps tmp/89o5k1229461042.png")
> system("convert tmp/92wq31229461042.ps tmp/92wq31229461042.png")
> system("convert tmp/10wg5q1229461043.ps tmp/10wg5q1229461043.png")
>
>
> proc.time()
user system elapsed
3.742 2.500 5.816