R version 2.8.0 (2008-10-20)
Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(94.6,116.1,95.9,107.5,104.7,116.7,102.8,112.5,98.1,113,113.9,126.4,80.9,114.1,95.7,112.5,113.2,112.4,105.9,113.1,108.8,116.3,102.3,111.7,99,118.8,100.7,116.5,115.5,125.1,100.7,113.1,109.9,119.6,114.6,114.4,85.4,114,100.5,117.8,114.8,117,116.5,120.9,112.9,115,102,117.3,106,119.4,105.3,114.9,118.8,125.8,106.1,117.6,109.3,117.6,117.2,114.9,92.5,121.9,104.2,117,112.5,106.4,122.4,110.5,113.3,113.6,100,114.2,110.7,125.4,112.8,124.6,109.8,120.2,117.3,120.8,109.1,111.4,115.9,124.1,96,120.2,99.8,125.5,116.8,116,115.7,117,99.4,105.7,94.3,102,91,106.4,93.2,96.9,103.1,107.6,94.1,98.8,91.8,101.1,102.7,105.7,82.6,104.6,89.1,103.2,104.5,101.6,105.1,106.7,95.1,99.5,88.7,101),dim=c(2,60),dimnames=list(c('T.I.P.','I.P.C.N.'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('T.I.P.','I.P.C.N.'),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 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '2'
> #'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
I.P.C.N. T.I.P. M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 116.1 94.6 1 0 0 0 0 0 0 0 0 0 0
2 107.5 95.9 0 1 0 0 0 0 0 0 0 0 0
3 116.7 104.7 0 0 1 0 0 0 0 0 0 0 0
4 112.5 102.8 0 0 0 1 0 0 0 0 0 0 0
5 113.0 98.1 0 0 0 0 1 0 0 0 0 0 0
6 126.4 113.9 0 0 0 0 0 1 0 0 0 0 0
7 114.1 80.9 0 0 0 0 0 0 1 0 0 0 0
8 112.5 95.7 0 0 0 0 0 0 0 1 0 0 0
9 112.4 113.2 0 0 0 0 0 0 0 0 1 0 0
10 113.1 105.9 0 0 0 0 0 0 0 0 0 1 0
11 116.3 108.8 0 0 0 0 0 0 0 0 0 0 1
12 111.7 102.3 0 0 0 0 0 0 0 0 0 0 0
13 118.8 99.0 1 0 0 0 0 0 0 0 0 0 0
14 116.5 100.7 0 1 0 0 0 0 0 0 0 0 0
15 125.1 115.5 0 0 1 0 0 0 0 0 0 0 0
16 113.1 100.7 0 0 0 1 0 0 0 0 0 0 0
17 119.6 109.9 0 0 0 0 1 0 0 0 0 0 0
18 114.4 114.6 0 0 0 0 0 1 0 0 0 0 0
19 114.0 85.4 0 0 0 0 0 0 1 0 0 0 0
20 117.8 100.5 0 0 0 0 0 0 0 1 0 0 0
21 117.0 114.8 0 0 0 0 0 0 0 0 1 0 0
22 120.9 116.5 0 0 0 0 0 0 0 0 0 1 0
23 115.0 112.9 0 0 0 0 0 0 0 0 0 0 1
24 117.3 102.0 0 0 0 0 0 0 0 0 0 0 0
25 119.4 106.0 1 0 0 0 0 0 0 0 0 0 0
26 114.9 105.3 0 1 0 0 0 0 0 0 0 0 0
27 125.8 118.8 0 0 1 0 0 0 0 0 0 0 0
28 117.6 106.1 0 0 0 1 0 0 0 0 0 0 0
29 117.6 109.3 0 0 0 0 1 0 0 0 0 0 0
30 114.9 117.2 0 0 0 0 0 1 0 0 0 0 0
31 121.9 92.5 0 0 0 0 0 0 1 0 0 0 0
32 117.0 104.2 0 0 0 0 0 0 0 1 0 0 0
33 106.4 112.5 0 0 0 0 0 0 0 0 1 0 0
34 110.5 122.4 0 0 0 0 0 0 0 0 0 1 0
35 113.6 113.3 0 0 0 0 0 0 0 0 0 0 1
36 114.2 100.0 0 0 0 0 0 0 0 0 0 0 0
37 125.4 110.7 1 0 0 0 0 0 0 0 0 0 0
38 124.6 112.8 0 1 0 0 0 0 0 0 0 0 0
39 120.2 109.8 0 0 1 0 0 0 0 0 0 0 0
40 120.8 117.3 0 0 0 1 0 0 0 0 0 0 0
41 111.4 109.1 0 0 0 0 1 0 0 0 0 0 0
42 124.1 115.9 0 0 0 0 0 1 0 0 0 0 0
43 120.2 96.0 0 0 0 0 0 0 1 0 0 0 0
44 125.5 99.8 0 0 0 0 0 0 0 1 0 0 0
45 116.0 116.8 0 0 0 0 0 0 0 0 1 0 0
46 117.0 115.7 0 0 0 0 0 0 0 0 0 1 0
47 105.7 99.4 0 0 0 0 0 0 0 0 0 0 1
48 102.0 94.3 0 0 0 0 0 0 0 0 0 0 0
49 106.4 91.0 1 0 0 0 0 0 0 0 0 0 0
50 96.9 93.2 0 1 0 0 0 0 0 0 0 0 0
51 107.6 103.1 0 0 1 0 0 0 0 0 0 0 0
52 98.8 94.1 0 0 0 1 0 0 0 0 0 0 0
53 101.1 91.8 0 0 0 0 1 0 0 0 0 0 0
54 105.7 102.7 0 0 0 0 0 1 0 0 0 0 0
55 104.6 82.6 0 0 0 0 0 0 1 0 0 0 0
56 103.2 89.1 0 0 0 0 0 0 0 1 0 0 0
57 101.6 104.5 0 0 0 0 0 0 0 0 1 0 0
58 106.7 105.1 0 0 0 0 0 0 0 0 0 1 0
59 99.5 95.1 0 0 0 0 0 0 0 0 0 0 1
60 101.0 88.7 0 0 0 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) T.I.P. M1 M2 M3 M4
23.6112 0.8786 5.5199 -0.7799 -1.5116 -2.6018
M5 M6 M7 M8 M9 M10
-2.1298 -5.6705 14.4885 5.6086 -11.6512 -9.3589
M11
-6.6354
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-11.2935 -2.7289 0.1934 2.6974 8.5955
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 23.61120 9.05787 2.607 0.012210 *
T.I.P. 0.87860 0.09067 9.690 8.76e-13 ***
M1 5.51991 2.82341 1.955 0.056539 .
M2 -0.77985 2.83668 -0.275 0.784585
M3 -1.51157 3.04624 -0.496 0.622060
M4 -2.60180 2.87762 -0.904 0.370527
M5 -2.12978 2.86726 -0.743 0.461305
M6 -5.67051 3.13959 -1.806 0.077303 .
M7 14.48847 2.95399 4.905 1.16e-05 ***
M8 5.60856 2.81221 1.994 0.051931 .
M9 -11.65121 3.11969 -3.735 0.000508 ***
M10 -9.35895 3.15014 -2.971 0.004666 **
M11 -6.63542 2.91425 -2.277 0.027388 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.446 on 47 degrees of freedom
Multiple R-squared: 0.7191, Adjusted R-squared: 0.6474
F-statistic: 10.03 on 12 and 47 DF, p-value: 2.472e-09
> 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.16074884 0.32149768 0.8392512
[2,] 0.09773836 0.19547673 0.9022616
[3,] 0.56841318 0.86317364 0.4315868
[4,] 0.46678731 0.93357462 0.5332127
[5,] 0.35078941 0.70157882 0.6492106
[6,] 0.29334112 0.58668223 0.7066589
[7,] 0.23828729 0.47657458 0.7617127
[8,] 0.19121159 0.38242319 0.8087884
[9,] 0.18957168 0.37914336 0.8104283
[10,] 0.15285559 0.30571117 0.8471444
[11,] 0.10176949 0.20353897 0.8982305
[12,] 0.06237260 0.12474519 0.9376274
[13,] 0.05471231 0.10942462 0.9452877
[14,] 0.03538157 0.07076314 0.9646184
[15,] 0.06786519 0.13573037 0.9321348
[16,] 0.06053698 0.12107396 0.9394630
[17,] 0.05467455 0.10934911 0.9453254
[18,] 0.07724949 0.15449898 0.9227505
[19,] 0.50944299 0.98111402 0.4905570
[20,] 0.54505753 0.90988495 0.4549425
[21,] 0.47257686 0.94515372 0.5274231
[22,] 0.40170227 0.80340453 0.5982977
[23,] 0.43663812 0.87327623 0.5633619
[24,] 0.40587557 0.81175115 0.5941244
[25,] 0.31626585 0.63253171 0.6837341
[26,] 0.61128714 0.77742573 0.3887129
[27,] 0.50938131 0.98123738 0.4906187
[28,] 0.37066561 0.74133121 0.6293344
[29,] 0.82206811 0.35586377 0.1779319
> postscript(file="/var/www/html/freestat/rcomp/tmp/19cbx1291212093.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/29cbx1291212093.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/39cbx1291212093.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/41lt11291212093.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/51lt11291212093.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
3.85290207 0.41047416 2.61047416 1.17004645 5.32746951 8.38625121
7 8 9 10 11 12
4.92121830 -0.80221405 0.98197213 5.80352525 3.73204664 -1.79244629
13 14 15 16 17 18
2.68704180 5.19317205 1.52154442 3.61511612 1.55993517 -4.22877201
19 20 21 22 23 24
0.86749758 0.28048384 4.17620476 4.29031643 -1.17023224 4.07113509
25 26 27 28 29 30
-2.86319044 -0.44840913 -0.67785078 3.37065125 0.08709793 -6.01314399
31 32 33 34 35 36
2.52940488 -3.77035320 -4.40300464 -11.29345074 -2.92167408 2.72834430
37 38 39 40 41 42
-0.99263208 2.66205633 1.62959067 -3.26972033 -5.93718115 4.32904200
43 44 45 46 47 48
-2.24571124 8.59550707 1.41899555 1.09320012 1.39092994 -4.46360945
49 50 51 52 53 54
-2.68412135 -7.81729341 -5.08375847 -4.88609348 -1.03732147 -2.47337721
55 56 57 58 59 60
-6.07240953 -4.30342366 -2.17416780 0.10640894 -1.03107026 -0.54342366
> postscript(file="/var/www/html/freestat/rcomp/tmp/61lt11291212093.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 3.85290207 NA
1 0.41047416 3.85290207
2 2.61047416 0.41047416
3 1.17004645 2.61047416
4 5.32746951 1.17004645
5 8.38625121 5.32746951
6 4.92121830 8.38625121
7 -0.80221405 4.92121830
8 0.98197213 -0.80221405
9 5.80352525 0.98197213
10 3.73204664 5.80352525
11 -1.79244629 3.73204664
12 2.68704180 -1.79244629
13 5.19317205 2.68704180
14 1.52154442 5.19317205
15 3.61511612 1.52154442
16 1.55993517 3.61511612
17 -4.22877201 1.55993517
18 0.86749758 -4.22877201
19 0.28048384 0.86749758
20 4.17620476 0.28048384
21 4.29031643 4.17620476
22 -1.17023224 4.29031643
23 4.07113509 -1.17023224
24 -2.86319044 4.07113509
25 -0.44840913 -2.86319044
26 -0.67785078 -0.44840913
27 3.37065125 -0.67785078
28 0.08709793 3.37065125
29 -6.01314399 0.08709793
30 2.52940488 -6.01314399
31 -3.77035320 2.52940488
32 -4.40300464 -3.77035320
33 -11.29345074 -4.40300464
34 -2.92167408 -11.29345074
35 2.72834430 -2.92167408
36 -0.99263208 2.72834430
37 2.66205633 -0.99263208
38 1.62959067 2.66205633
39 -3.26972033 1.62959067
40 -5.93718115 -3.26972033
41 4.32904200 -5.93718115
42 -2.24571124 4.32904200
43 8.59550707 -2.24571124
44 1.41899555 8.59550707
45 1.09320012 1.41899555
46 1.39092994 1.09320012
47 -4.46360945 1.39092994
48 -2.68412135 -4.46360945
49 -7.81729341 -2.68412135
50 -5.08375847 -7.81729341
51 -4.88609348 -5.08375847
52 -1.03732147 -4.88609348
53 -2.47337721 -1.03732147
54 -6.07240953 -2.47337721
55 -4.30342366 -6.07240953
56 -2.17416780 -4.30342366
57 0.10640894 -2.17416780
58 -1.03107026 0.10640894
59 -0.54342366 -1.03107026
60 NA -0.54342366
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.41047416 3.85290207
[2,] 2.61047416 0.41047416
[3,] 1.17004645 2.61047416
[4,] 5.32746951 1.17004645
[5,] 8.38625121 5.32746951
[6,] 4.92121830 8.38625121
[7,] -0.80221405 4.92121830
[8,] 0.98197213 -0.80221405
[9,] 5.80352525 0.98197213
[10,] 3.73204664 5.80352525
[11,] -1.79244629 3.73204664
[12,] 2.68704180 -1.79244629
[13,] 5.19317205 2.68704180
[14,] 1.52154442 5.19317205
[15,] 3.61511612 1.52154442
[16,] 1.55993517 3.61511612
[17,] -4.22877201 1.55993517
[18,] 0.86749758 -4.22877201
[19,] 0.28048384 0.86749758
[20,] 4.17620476 0.28048384
[21,] 4.29031643 4.17620476
[22,] -1.17023224 4.29031643
[23,] 4.07113509 -1.17023224
[24,] -2.86319044 4.07113509
[25,] -0.44840913 -2.86319044
[26,] -0.67785078 -0.44840913
[27,] 3.37065125 -0.67785078
[28,] 0.08709793 3.37065125
[29,] -6.01314399 0.08709793
[30,] 2.52940488 -6.01314399
[31,] -3.77035320 2.52940488
[32,] -4.40300464 -3.77035320
[33,] -11.29345074 -4.40300464
[34,] -2.92167408 -11.29345074
[35,] 2.72834430 -2.92167408
[36,] -0.99263208 2.72834430
[37,] 2.66205633 -0.99263208
[38,] 1.62959067 2.66205633
[39,] -3.26972033 1.62959067
[40,] -5.93718115 -3.26972033
[41,] 4.32904200 -5.93718115
[42,] -2.24571124 4.32904200
[43,] 8.59550707 -2.24571124
[44,] 1.41899555 8.59550707
[45,] 1.09320012 1.41899555
[46,] 1.39092994 1.09320012
[47,] -4.46360945 1.39092994
[48,] -2.68412135 -4.46360945
[49,] -7.81729341 -2.68412135
[50,] -5.08375847 -7.81729341
[51,] -4.88609348 -5.08375847
[52,] -1.03732147 -4.88609348
[53,] -2.47337721 -1.03732147
[54,] -6.07240953 -2.47337721
[55,] -4.30342366 -6.07240953
[56,] -2.17416780 -4.30342366
[57,] 0.10640894 -2.17416780
[58,] -1.03107026 0.10640894
[59,] -0.54342366 -1.03107026
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.41047416 3.85290207
2 2.61047416 0.41047416
3 1.17004645 2.61047416
4 5.32746951 1.17004645
5 8.38625121 5.32746951
6 4.92121830 8.38625121
7 -0.80221405 4.92121830
8 0.98197213 -0.80221405
9 5.80352525 0.98197213
10 3.73204664 5.80352525
11 -1.79244629 3.73204664
12 2.68704180 -1.79244629
13 5.19317205 2.68704180
14 1.52154442 5.19317205
15 3.61511612 1.52154442
16 1.55993517 3.61511612
17 -4.22877201 1.55993517
18 0.86749758 -4.22877201
19 0.28048384 0.86749758
20 4.17620476 0.28048384
21 4.29031643 4.17620476
22 -1.17023224 4.29031643
23 4.07113509 -1.17023224
24 -2.86319044 4.07113509
25 -0.44840913 -2.86319044
26 -0.67785078 -0.44840913
27 3.37065125 -0.67785078
28 0.08709793 3.37065125
29 -6.01314399 0.08709793
30 2.52940488 -6.01314399
31 -3.77035320 2.52940488
32 -4.40300464 -3.77035320
33 -11.29345074 -4.40300464
34 -2.92167408 -11.29345074
35 2.72834430 -2.92167408
36 -0.99263208 2.72834430
37 2.66205633 -0.99263208
38 1.62959067 2.66205633
39 -3.26972033 1.62959067
40 -5.93718115 -3.26972033
41 4.32904200 -5.93718115
42 -2.24571124 4.32904200
43 8.59550707 -2.24571124
44 1.41899555 8.59550707
45 1.09320012 1.41899555
46 1.39092994 1.09320012
47 -4.46360945 1.39092994
48 -2.68412135 -4.46360945
49 -7.81729341 -2.68412135
50 -5.08375847 -7.81729341
51 -4.88609348 -5.08375847
52 -1.03732147 -4.88609348
53 -2.47337721 -1.03732147
54 -6.07240953 -2.47337721
55 -4.30342366 -6.07240953
56 -2.17416780 -4.30342366
57 0.10640894 -2.17416780
58 -1.03107026 0.10640894
59 -0.54342366 -1.03107026
> 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/7cual1291212093.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/853ro1291212093.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/953ro1291212093.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/1053ro1291212093.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/freestat/rcomp/tmp/1184pu1291212093.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/12tmo01291212093.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/13in3t1291212093.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/14tfkf1291212093.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/1576051291212093.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/16apyt1291212093.tab")
+ }
>
> try(system("convert tmp/19cbx1291212093.ps tmp/19cbx1291212093.png",intern=TRUE))
character(0)
> try(system("convert tmp/29cbx1291212093.ps tmp/29cbx1291212093.png",intern=TRUE))
character(0)
> try(system("convert tmp/39cbx1291212093.ps tmp/39cbx1291212093.png",intern=TRUE))
character(0)
> try(system("convert tmp/41lt11291212093.ps tmp/41lt11291212093.png",intern=TRUE))
character(0)
> try(system("convert tmp/51lt11291212093.ps tmp/51lt11291212093.png",intern=TRUE))
character(0)
> try(system("convert tmp/61lt11291212093.ps tmp/61lt11291212093.png",intern=TRUE))
character(0)
> try(system("convert tmp/7cual1291212093.ps tmp/7cual1291212093.png",intern=TRUE))
character(0)
> try(system("convert tmp/853ro1291212093.ps tmp/853ro1291212093.png",intern=TRUE))
character(0)
> try(system("convert tmp/953ro1291212093.ps tmp/953ro1291212093.png",intern=TRUE))
character(0)
> try(system("convert tmp/1053ro1291212093.ps tmp/1053ro1291212093.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.768 2.475 4.084