R version 2.12.0 (2010-10-15)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-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(112.3,1,117.3,1,111.1,1,102.2,1,104.3,1,122.9,0,107.6,0,121.3,0,131.5,0,89,0,104.4,0,128.9,0,135.9,0,133.3,0,121.3,0,120.5,0,120.4,0,137.9,0,126.1,0,133.2,0,151.1,0,105,0,119,0,140.4,0,156.6,1,137.1,1,122.7,1,125.8,1,139.3,1,134.9,1,149.2,1,132.3,1,149,1,117.2,1,119.6,1,152,1,149.4,1,127.3,1,114.1,1,102.1,1,107.7,1,104.4,1,102.1,1,96,1,109.3,1,90,1,83.9,1,112,1,114.3,1,103.6,1,91.7,1,80.8,1,87.2,1,109.2,1,102.7,1,95.1,1,117.5,1,85.1,1,92.1,1,113.5,1),dim=c(2,60),dimnames=list(c('Promet','Dummy'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Promet','Dummy'),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
Promet Dummy M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 112.3 1 1 0 0 0 0 0 0 0 0 0 0 1
2 117.3 1 0 1 0 0 0 0 0 0 0 0 0 2
3 111.1 1 0 0 1 0 0 0 0 0 0 0 0 3
4 102.2 1 0 0 0 1 0 0 0 0 0 0 0 4
5 104.3 1 0 0 0 0 1 0 0 0 0 0 0 5
6 122.9 0 0 0 0 0 0 1 0 0 0 0 0 6
7 107.6 0 0 0 0 0 0 0 1 0 0 0 0 7
8 121.3 0 0 0 0 0 0 0 0 1 0 0 0 8
9 131.5 0 0 0 0 0 0 0 0 0 1 0 0 9
10 89.0 0 0 0 0 0 0 0 0 0 0 1 0 10
11 104.4 0 0 0 0 0 0 0 0 0 0 0 1 11
12 128.9 0 0 0 0 0 0 0 0 0 0 0 0 12
13 135.9 0 1 0 0 0 0 0 0 0 0 0 0 13
14 133.3 0 0 1 0 0 0 0 0 0 0 0 0 14
15 121.3 0 0 0 1 0 0 0 0 0 0 0 0 15
16 120.5 0 0 0 0 1 0 0 0 0 0 0 0 16
17 120.4 0 0 0 0 0 1 0 0 0 0 0 0 17
18 137.9 0 0 0 0 0 0 1 0 0 0 0 0 18
19 126.1 0 0 0 0 0 0 0 1 0 0 0 0 19
20 133.2 0 0 0 0 0 0 0 0 1 0 0 0 20
21 151.1 0 0 0 0 0 0 0 0 0 1 0 0 21
22 105.0 0 0 0 0 0 0 0 0 0 0 1 0 22
23 119.0 0 0 0 0 0 0 0 0 0 0 0 1 23
24 140.4 0 0 0 0 0 0 0 0 0 0 0 0 24
25 156.6 1 1 0 0 0 0 0 0 0 0 0 0 25
26 137.1 1 0 1 0 0 0 0 0 0 0 0 0 26
27 122.7 1 0 0 1 0 0 0 0 0 0 0 0 27
28 125.8 1 0 0 0 1 0 0 0 0 0 0 0 28
29 139.3 1 0 0 0 0 1 0 0 0 0 0 0 29
30 134.9 1 0 0 0 0 0 1 0 0 0 0 0 30
31 149.2 1 0 0 0 0 0 0 1 0 0 0 0 31
32 132.3 1 0 0 0 0 0 0 0 1 0 0 0 32
33 149.0 1 0 0 0 0 0 0 0 0 1 0 0 33
34 117.2 1 0 0 0 0 0 0 0 0 0 1 0 34
35 119.6 1 0 0 0 0 0 0 0 0 0 0 1 35
36 152.0 1 0 0 0 0 0 0 0 0 0 0 0 36
37 149.4 1 1 0 0 0 0 0 0 0 0 0 0 37
38 127.3 1 0 1 0 0 0 0 0 0 0 0 0 38
39 114.1 1 0 0 1 0 0 0 0 0 0 0 0 39
40 102.1 1 0 0 0 1 0 0 0 0 0 0 0 40
41 107.7 1 0 0 0 0 1 0 0 0 0 0 0 41
42 104.4 1 0 0 0 0 0 1 0 0 0 0 0 42
43 102.1 1 0 0 0 0 0 0 1 0 0 0 0 43
44 96.0 1 0 0 0 0 0 0 0 1 0 0 0 44
45 109.3 1 0 0 0 0 0 0 0 0 1 0 0 45
46 90.0 1 0 0 0 0 0 0 0 0 0 1 0 46
47 83.9 1 0 0 0 0 0 0 0 0 0 0 1 47
48 112.0 1 0 0 0 0 0 0 0 0 0 0 0 48
49 114.3 1 1 0 0 0 0 0 0 0 0 0 0 49
50 103.6 1 0 1 0 0 0 0 0 0 0 0 0 50
51 91.7 1 0 0 1 0 0 0 0 0 0 0 0 51
52 80.8 1 0 0 0 1 0 0 0 0 0 0 0 52
53 87.2 1 0 0 0 0 1 0 0 0 0 0 0 53
54 109.2 1 0 0 0 0 0 1 0 0 0 0 0 54
55 102.7 1 0 0 0 0 0 0 1 0 0 0 0 55
56 95.1 1 0 0 0 0 0 0 0 1 0 0 0 56
57 117.5 1 0 0 0 0 0 0 0 0 1 0 0 57
58 85.1 1 0 0 0 0 0 0 0 0 0 1 0 58
59 92.1 1 0 0 0 0 0 0 0 0 0 0 1 59
60 113.5 1 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) Dummy M1 M2 M3 M4
143.4659 -0.7818 0.3295 -9.2717 -20.4329 -25.9541
M5 M6 M7 M8 M9 M10
-20.0753 -9.7728 -13.7140 -15.2952 1.1836 -32.8576
M11 t
-25.9388 -0.3788
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-30.335 -10.923 -2.579 10.764 31.973
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 143.4659 8.2319 17.428 < 2e-16 ***
Dummy -0.7818 6.0306 -0.130 0.89742
M1 0.3295 10.2925 0.032 0.97460
M2 -9.2717 10.2531 -0.904 0.37056
M3 -20.4329 10.2160 -2.000 0.05142 .
M4 -25.9541 10.1814 -2.549 0.01420 *
M5 -20.0753 10.1492 -1.978 0.05394 .
M6 -9.7728 9.9687 -0.980 0.33204
M7 -13.7140 9.9543 -1.378 0.17496
M8 -15.2952 9.9425 -1.538 0.13081
M9 1.1836 9.9333 0.119 0.90567
M10 -32.8576 9.9267 -3.310 0.00182 **
M11 -25.9388 9.9228 -2.614 0.01205 *
t -0.3788 0.1615 -2.345 0.02341 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15.69 on 46 degrees of freedom
Multiple R-squared: 0.4668, Adjusted R-squared: 0.3161
F-statistic: 3.098 on 13 and 46 DF, p-value: 0.002306
> if (n > n25) {
+ kp3 <- k + 3
+ nmkm3 <- n - k - 3
+ gqarr <- array(NA, dim=c(nmkm3-kp3+1,3))
+ numgqtests <- 0
+ numsignificant1 <- 0
+ numsignificant5 <- 0
+ numsignificant10 <- 0
+ for (mypoint in kp3:nmkm3) {
+ j <- 0
+ numgqtests <- numgqtests + 1
+ for (myalt in c('greater', 'two.sided', 'less')) {
+ j <- j + 1
+ gqarr[mypoint-kp3+1,j] <- gqtest(mylm, point=mypoint, alternative=myalt)$p.value
+ }
+ if (gqarr[mypoint-kp3+1,2] < 0.01) numsignificant1 <- numsignificant1 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.05) numsignificant5 <- numsignificant5 + 1
+ if (gqarr[mypoint-kp3+1,2] < 0.10) numsignificant10 <- numsignificant10 + 1
+ }
+ gqarr
+ }
[,1] [,2] [,3]
[1,] 8.143589e-02 0.1628717872 0.91856411
[2,] 2.605456e-02 0.0521091156 0.97394544
[3,] 9.103466e-03 0.0182069319 0.99089653
[4,] 3.831748e-03 0.0076634967 0.99616825
[5,] 1.663717e-03 0.0033274330 0.99833628
[6,] 4.743336e-04 0.0009486671 0.99952567
[7,] 1.286189e-04 0.0002572378 0.99987138
[8,] 5.230331e-05 0.0001046066 0.99994770
[9,] 1.002242e-04 0.0002004483 0.99989978
[10,] 8.187188e-04 0.0016374377 0.99918128
[11,] 3.187167e-03 0.0063743344 0.99681283
[12,] 1.554452e-03 0.0031089043 0.99844555
[13,] 1.698825e-03 0.0033976493 0.99830118
[14,] 2.065318e-03 0.0041306355 0.99793468
[15,] 9.781531e-03 0.0195630618 0.99021847
[16,] 1.225501e-02 0.0245100238 0.98774499
[17,] 1.079591e-02 0.0215918265 0.98920409
[18,] 6.938788e-03 0.0138775768 0.99306121
[19,] 5.724158e-03 0.0114483164 0.99427584
[20,] 1.152803e-02 0.0230560539 0.98847197
[21,] 4.853435e-02 0.0970686965 0.95146565
[22,] 2.131932e-01 0.4263863136 0.78680684
[23,] 4.582713e-01 0.9165426537 0.54172867
[24,] 7.872423e-01 0.4255153877 0.21275769
[25,] 9.861316e-01 0.0277367715 0.01386839
[26,] 9.773070e-01 0.0453860355 0.02269302
[27,] 9.440961e-01 0.1118078516 0.05590393
> postscript(file="/var/www/rcomp/tmp/1net51292963184.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/rcomp/tmp/2yob81292963184.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/rcomp/tmp/3yob81292963184.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/rcomp/tmp/4yob81292963184.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/rcomp/tmp/5yob81292963184.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
-30.3348768 -15.3548768 -10.0148768 -13.0148768 -16.4148768 -8.5202956
7 8 9 10 11 12
-19.5002956 -3.8402956 -9.7402956 -17.8202956 -8.9602956 -10.0202956
13 14 15 16 17 18
-2.9710345 4.4089655 3.9489655 9.0489655 3.4489655 11.0253202
19 20 21 22 23 24
3.5453202 12.6053202 14.4053202 2.7253202 10.1853202 6.0253202
25 26 27 28 29 30
23.0563547 13.5363547 10.6763547 19.6763547 27.6763547 13.3527094
31 32 33 34 35 36
31.9727094 17.0327094 17.6327094 20.2527094 16.1127094 22.9527094
37 38 39 40 41 42
20.4019704 8.2819704 6.6219704 0.5219704 0.6219704 -12.6016749
43 44 45 46 47 48
-10.5816749 -14.7216749 -17.5216749 -2.4016749 -15.0416749 -12.5016749
49 50 51 52 53 54
-10.1524138 -10.8724138 -11.2324138 -16.2324138 -15.3324138 -3.2560591
55 56 57 58 59 60
-5.4360591 -11.0760591 -4.7760591 -2.7560591 -2.2960591 -6.4560591
> postscript(file="/var/www/rcomp/tmp/6qfsa1292963184.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 -30.3348768 NA
1 -15.3548768 -30.3348768
2 -10.0148768 -15.3548768
3 -13.0148768 -10.0148768
4 -16.4148768 -13.0148768
5 -8.5202956 -16.4148768
6 -19.5002956 -8.5202956
7 -3.8402956 -19.5002956
8 -9.7402956 -3.8402956
9 -17.8202956 -9.7402956
10 -8.9602956 -17.8202956
11 -10.0202956 -8.9602956
12 -2.9710345 -10.0202956
13 4.4089655 -2.9710345
14 3.9489655 4.4089655
15 9.0489655 3.9489655
16 3.4489655 9.0489655
17 11.0253202 3.4489655
18 3.5453202 11.0253202
19 12.6053202 3.5453202
20 14.4053202 12.6053202
21 2.7253202 14.4053202
22 10.1853202 2.7253202
23 6.0253202 10.1853202
24 23.0563547 6.0253202
25 13.5363547 23.0563547
26 10.6763547 13.5363547
27 19.6763547 10.6763547
28 27.6763547 19.6763547
29 13.3527094 27.6763547
30 31.9727094 13.3527094
31 17.0327094 31.9727094
32 17.6327094 17.0327094
33 20.2527094 17.6327094
34 16.1127094 20.2527094
35 22.9527094 16.1127094
36 20.4019704 22.9527094
37 8.2819704 20.4019704
38 6.6219704 8.2819704
39 0.5219704 6.6219704
40 0.6219704 0.5219704
41 -12.6016749 0.6219704
42 -10.5816749 -12.6016749
43 -14.7216749 -10.5816749
44 -17.5216749 -14.7216749
45 -2.4016749 -17.5216749
46 -15.0416749 -2.4016749
47 -12.5016749 -15.0416749
48 -10.1524138 -12.5016749
49 -10.8724138 -10.1524138
50 -11.2324138 -10.8724138
51 -16.2324138 -11.2324138
52 -15.3324138 -16.2324138
53 -3.2560591 -15.3324138
54 -5.4360591 -3.2560591
55 -11.0760591 -5.4360591
56 -4.7760591 -11.0760591
57 -2.7560591 -4.7760591
58 -2.2960591 -2.7560591
59 -6.4560591 -2.2960591
60 NA -6.4560591
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -15.3548768 -30.3348768
[2,] -10.0148768 -15.3548768
[3,] -13.0148768 -10.0148768
[4,] -16.4148768 -13.0148768
[5,] -8.5202956 -16.4148768
[6,] -19.5002956 -8.5202956
[7,] -3.8402956 -19.5002956
[8,] -9.7402956 -3.8402956
[9,] -17.8202956 -9.7402956
[10,] -8.9602956 -17.8202956
[11,] -10.0202956 -8.9602956
[12,] -2.9710345 -10.0202956
[13,] 4.4089655 -2.9710345
[14,] 3.9489655 4.4089655
[15,] 9.0489655 3.9489655
[16,] 3.4489655 9.0489655
[17,] 11.0253202 3.4489655
[18,] 3.5453202 11.0253202
[19,] 12.6053202 3.5453202
[20,] 14.4053202 12.6053202
[21,] 2.7253202 14.4053202
[22,] 10.1853202 2.7253202
[23,] 6.0253202 10.1853202
[24,] 23.0563547 6.0253202
[25,] 13.5363547 23.0563547
[26,] 10.6763547 13.5363547
[27,] 19.6763547 10.6763547
[28,] 27.6763547 19.6763547
[29,] 13.3527094 27.6763547
[30,] 31.9727094 13.3527094
[31,] 17.0327094 31.9727094
[32,] 17.6327094 17.0327094
[33,] 20.2527094 17.6327094
[34,] 16.1127094 20.2527094
[35,] 22.9527094 16.1127094
[36,] 20.4019704 22.9527094
[37,] 8.2819704 20.4019704
[38,] 6.6219704 8.2819704
[39,] 0.5219704 6.6219704
[40,] 0.6219704 0.5219704
[41,] -12.6016749 0.6219704
[42,] -10.5816749 -12.6016749
[43,] -14.7216749 -10.5816749
[44,] -17.5216749 -14.7216749
[45,] -2.4016749 -17.5216749
[46,] -15.0416749 -2.4016749
[47,] -12.5016749 -15.0416749
[48,] -10.1524138 -12.5016749
[49,] -10.8724138 -10.1524138
[50,] -11.2324138 -10.8724138
[51,] -16.2324138 -11.2324138
[52,] -15.3324138 -16.2324138
[53,] -3.2560591 -15.3324138
[54,] -5.4360591 -3.2560591
[55,] -11.0760591 -5.4360591
[56,] -4.7760591 -11.0760591
[57,] -2.7560591 -4.7760591
[58,] -2.2960591 -2.7560591
[59,] -6.4560591 -2.2960591
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -15.3548768 -30.3348768
2 -10.0148768 -15.3548768
3 -13.0148768 -10.0148768
4 -16.4148768 -13.0148768
5 -8.5202956 -16.4148768
6 -19.5002956 -8.5202956
7 -3.8402956 -19.5002956
8 -9.7402956 -3.8402956
9 -17.8202956 -9.7402956
10 -8.9602956 -17.8202956
11 -10.0202956 -8.9602956
12 -2.9710345 -10.0202956
13 4.4089655 -2.9710345
14 3.9489655 4.4089655
15 9.0489655 3.9489655
16 3.4489655 9.0489655
17 11.0253202 3.4489655
18 3.5453202 11.0253202
19 12.6053202 3.5453202
20 14.4053202 12.6053202
21 2.7253202 14.4053202
22 10.1853202 2.7253202
23 6.0253202 10.1853202
24 23.0563547 6.0253202
25 13.5363547 23.0563547
26 10.6763547 13.5363547
27 19.6763547 10.6763547
28 27.6763547 19.6763547
29 13.3527094 27.6763547
30 31.9727094 13.3527094
31 17.0327094 31.9727094
32 17.6327094 17.0327094
33 20.2527094 17.6327094
34 16.1127094 20.2527094
35 22.9527094 16.1127094
36 20.4019704 22.9527094
37 8.2819704 20.4019704
38 6.6219704 8.2819704
39 0.5219704 6.6219704
40 0.6219704 0.5219704
41 -12.6016749 0.6219704
42 -10.5816749 -12.6016749
43 -14.7216749 -10.5816749
44 -17.5216749 -14.7216749
45 -2.4016749 -17.5216749
46 -15.0416749 -2.4016749
47 -12.5016749 -15.0416749
48 -10.1524138 -12.5016749
49 -10.8724138 -10.1524138
50 -11.2324138 -10.8724138
51 -16.2324138 -11.2324138
52 -15.3324138 -16.2324138
53 -3.2560591 -15.3324138
54 -5.4360591 -3.2560591
55 -11.0760591 -5.4360591
56 -4.7760591 -11.0760591
57 -2.7560591 -4.7760591
58 -2.2960591 -2.7560591
59 -6.4560591 -2.2960591
> 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/rcomp/tmp/7169v1292963184.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/rcomp/tmp/8169v1292963184.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/rcomp/tmp/9169v1292963184.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/rcomp/tmp/10uxqg1292963184.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/11xg741292963184.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/rcomp/tmp/12iy5a1292963184.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/rcomp/tmp/13eql11292963184.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/rcomp/tmp/14i9271292963184.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/rcomp/tmp/15l9iv1292963184.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/rcomp/tmp/167azi1292963184.tab")
+ }
>
> try(system("convert tmp/1net51292963184.ps tmp/1net51292963184.png",intern=TRUE))
character(0)
> try(system("convert tmp/2yob81292963184.ps tmp/2yob81292963184.png",intern=TRUE))
character(0)
> try(system("convert tmp/3yob81292963184.ps tmp/3yob81292963184.png",intern=TRUE))
character(0)
> try(system("convert tmp/4yob81292963184.ps tmp/4yob81292963184.png",intern=TRUE))
character(0)
> try(system("convert tmp/5yob81292963184.ps tmp/5yob81292963184.png",intern=TRUE))
character(0)
> try(system("convert tmp/6qfsa1292963184.ps tmp/6qfsa1292963184.png",intern=TRUE))
character(0)
> try(system("convert tmp/7169v1292963184.ps tmp/7169v1292963184.png",intern=TRUE))
character(0)
> try(system("convert tmp/8169v1292963184.ps tmp/8169v1292963184.png",intern=TRUE))
character(0)
> try(system("convert tmp/9169v1292963184.ps tmp/9169v1292963184.png",intern=TRUE))
character(0)
> try(system("convert tmp/10uxqg1292963184.ps tmp/10uxqg1292963184.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.080 1.330 4.395