R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(110.5,55,110.8,48.7,104.2,70.3,88.9,94.8,89.8,58.5,90,62.4,93.9,56.7,91.3,65.1,87.8,114.4,99.7,50.7,73.5,44.5,79.2,72,96.9,61.2,95.2,68.4,95.6,78.7,89.7,64.1,92.8,64.6,88,71.9,101.1,71,92.7,76.4,95.8,117.3,103.8,66.1,81.8,57.3,87.1,75,105.9,63.8,108.1,62.2,102.6,75.4,93.7,58,103.5,62.1,100.6,99.2,113.3,70.7,102.4,73.3,102.1,111.2,106.9,68.9,87.3,57.6,93.1,72.9,109.1,75.9,120.3,79.4,104.9,96.9,92.6,75.2,109.8,60.3,111.4,88.9,117.9,90.5,121.6,79.9,117.8,116.3,124.2,95.2,106.8,81.5,102.7,89.1,116.8,76,113.6,100.5,96.1,83.9,85,75.1,83.2,69.5,84.9,95.1,83,90.1,79.6,78.4,83.2,113.8,83.8,73.6,82.8,56.5,71.4,97.7),dim=c(2,60),dimnames=list(c('prod','inv
'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('prod','inv
'),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 = '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
prod inv\r M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 110.5 55.0 1 0 0 0 0 0 0 0 0 0 0
2 110.8 48.7 0 1 0 0 0 0 0 0 0 0 0
3 104.2 70.3 0 0 1 0 0 0 0 0 0 0 0
4 88.9 94.8 0 0 0 1 0 0 0 0 0 0 0
5 89.8 58.5 0 0 0 0 1 0 0 0 0 0 0
6 90.0 62.4 0 0 0 0 0 1 0 0 0 0 0
7 93.9 56.7 0 0 0 0 0 0 1 0 0 0 0
8 91.3 65.1 0 0 0 0 0 0 0 1 0 0 0
9 87.8 114.4 0 0 0 0 0 0 0 0 1 0 0
10 99.7 50.7 0 0 0 0 0 0 0 0 0 1 0
11 73.5 44.5 0 0 0 0 0 0 0 0 0 0 1
12 79.2 72.0 0 0 0 0 0 0 0 0 0 0 0
13 96.9 61.2 1 0 0 0 0 0 0 0 0 0 0
14 95.2 68.4 0 1 0 0 0 0 0 0 0 0 0
15 95.6 78.7 0 0 1 0 0 0 0 0 0 0 0
16 89.7 64.1 0 0 0 1 0 0 0 0 0 0 0
17 92.8 64.6 0 0 0 0 1 0 0 0 0 0 0
18 88.0 71.9 0 0 0 0 0 1 0 0 0 0 0
19 101.1 71.0 0 0 0 0 0 0 1 0 0 0 0
20 92.7 76.4 0 0 0 0 0 0 0 1 0 0 0
21 95.8 117.3 0 0 0 0 0 0 0 0 1 0 0
22 103.8 66.1 0 0 0 0 0 0 0 0 0 1 0
23 81.8 57.3 0 0 0 0 0 0 0 0 0 0 1
24 87.1 75.0 0 0 0 0 0 0 0 0 0 0 0
25 105.9 63.8 1 0 0 0 0 0 0 0 0 0 0
26 108.1 62.2 0 1 0 0 0 0 0 0 0 0 0
27 102.6 75.4 0 0 1 0 0 0 0 0 0 0 0
28 93.7 58.0 0 0 0 1 0 0 0 0 0 0 0
29 103.5 62.1 0 0 0 0 1 0 0 0 0 0 0
30 100.6 99.2 0 0 0 0 0 1 0 0 0 0 0
31 113.3 70.7 0 0 0 0 0 0 1 0 0 0 0
32 102.4 73.3 0 0 0 0 0 0 0 1 0 0 0
33 102.1 111.2 0 0 0 0 0 0 0 0 1 0 0
34 106.9 68.9 0 0 0 0 0 0 0 0 0 1 0
35 87.3 57.6 0 0 0 0 0 0 0 0 0 0 1
36 93.1 72.9 0 0 0 0 0 0 0 0 0 0 0
37 109.1 75.9 1 0 0 0 0 0 0 0 0 0 0
38 120.3 79.4 0 1 0 0 0 0 0 0 0 0 0
39 104.9 96.9 0 0 1 0 0 0 0 0 0 0 0
40 92.6 75.2 0 0 0 1 0 0 0 0 0 0 0
41 109.8 60.3 0 0 0 0 1 0 0 0 0 0 0
42 111.4 88.9 0 0 0 0 0 1 0 0 0 0 0
43 117.9 90.5 0 0 0 0 0 0 1 0 0 0 0
44 121.6 79.9 0 0 0 0 0 0 0 1 0 0 0
45 117.8 116.3 0 0 0 0 0 0 0 0 1 0 0
46 124.2 95.2 0 0 0 0 0 0 0 0 0 1 0
47 106.8 81.5 0 0 0 0 0 0 0 0 0 0 1
48 102.7 89.1 0 0 0 0 0 0 0 0 0 0 0
49 116.8 76.0 1 0 0 0 0 0 0 0 0 0 0
50 113.6 100.5 0 1 0 0 0 0 0 0 0 0 0
51 96.1 83.9 0 0 1 0 0 0 0 0 0 0 0
52 85.0 75.1 0 0 0 1 0 0 0 0 0 0 0
53 83.2 69.5 0 0 0 0 1 0 0 0 0 0 0
54 84.9 95.1 0 0 0 0 0 1 0 0 0 0 0
55 83.0 90.1 0 0 0 0 0 0 1 0 0 0 0
56 79.6 78.4 0 0 0 0 0 0 0 1 0 0 0
57 83.2 113.8 0 0 0 0 0 0 0 0 1 0 0
58 83.8 73.6 0 0 0 0 0 0 0 0 0 1 0
59 82.8 56.5 0 0 0 0 0 0 0 0 0 0 1
60 71.4 97.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) `inv\r` M1 M2 M3 M4
68.0313 0.2295 24.5735 25.0804 14.0489 5.0932
M5 M6 M7 M8 M9 M10
13.3293 7.7842 16.4115 12.3623 3.0064 19.3761
M11
4.7572
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-22.1221 -5.2305 0.4708 6.1165 22.8682
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 68.0313 11.6171 5.856 4.44e-07 ***
`inv\r` 0.2295 0.1291 1.778 0.081896 .
M1 24.5735 7.2882 3.372 0.001502 **
M2 25.0804 7.1338 3.516 0.000983 ***
M3 14.0489 7.0277 1.999 0.051403 .
M4 5.0932 7.1012 0.717 0.476786
M5 13.3293 7.4157 1.797 0.078693 .
M6 7.7842 7.0331 1.107 0.274015
M7 16.4115 7.0639 2.323 0.024541 *
M8 12.3623 7.0810 1.746 0.087370 .
M9 3.0064 8.2355 0.365 0.716713
M10 19.3761 7.1557 2.708 0.009415 **
M11 4.7572 7.5731 0.628 0.532932
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11.11 on 47 degrees of freedom
Multiple R-squared: 0.3771, Adjusted R-squared: 0.2181
F-statistic: 2.372 on 12 and 47 DF, p-value: 0.01753
> 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,] 2.098870e-01 4.197740e-01 0.7901130
[2,] 1.073789e-01 2.147578e-01 0.8926211
[3,] 4.617049e-02 9.234098e-02 0.9538295
[4,] 3.634500e-02 7.269000e-02 0.9636550
[5,] 1.652181e-02 3.304361e-02 0.9834782
[6,] 1.013251e-02 2.026501e-02 0.9898675
[7,] 5.052797e-03 1.010559e-02 0.9949472
[8,] 3.525815e-03 7.051630e-03 0.9964742
[9,] 1.965054e-03 3.930108e-03 0.9980349
[10,] 7.771923e-04 1.554385e-03 0.9992228
[11,] 3.576332e-04 7.152663e-04 0.9996424
[12,] 1.302792e-04 2.605584e-04 0.9998697
[13,] 4.757304e-05 9.514609e-05 0.9999524
[14,] 6.293260e-05 1.258652e-04 0.9999371
[15,] 6.053531e-05 1.210706e-04 0.9999395
[16,] 1.296857e-04 2.593715e-04 0.9998703
[17,] 8.560865e-05 1.712173e-04 0.9999144
[18,] 5.775432e-05 1.155086e-04 0.9999422
[19,] 2.320627e-05 4.641254e-05 0.9999768
[20,] 1.209723e-05 2.419447e-05 0.9999879
[21,] 8.614948e-06 1.722990e-05 0.9999914
[22,] 2.915315e-06 5.830631e-06 0.9999971
[23,] 3.274506e-06 6.549011e-06 0.9999967
[24,] 9.216113e-07 1.843223e-06 0.9999991
[25,] 2.569257e-07 5.138515e-07 0.9999997
[26,] 1.107920e-06 2.215841e-06 0.9999989
[27,] 5.279908e-06 1.055982e-05 0.9999947
[28,] 1.434179e-05 2.868358e-05 0.9999857
[29,] 5.738561e-04 1.147712e-03 0.9994261
> postscript(file="/var/www/html/rcomp/tmp/19zj41258637484.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/2gkm41258637484.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/3ag0g1258637484.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/4cp4j1258637484.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/50xb01258637484.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 = 60
Frequency = 1
1 2 3 4 5 6
5.2718739 6.5109633 5.9849847 -5.9824277 -4.9871852 -0.1372460
7 8 9 10 11 12
-3.5562749 -4.0350229 -9.4940971 0.6561910 -9.5018742 -5.3563355
13 14 15 16 17 18
-9.7511154 -13.6104704 -4.5429363 1.8636645 -3.3872230 -4.3176329
19 20 21 22 23 24
0.3616691 -5.2285356 -2.1596889 1.2216691 -4.1396586 1.8551213
25 26 27 28 29 30
-1.3478528 0.7125189 3.2144612 7.2637024 7.8865630 2.0166238
31 32 33 34 35 36
12.6305235 5.1829590 5.5403490 3.6790288 1.2914871 8.3371016
37 38 39 40 41 42
-0.9249771 8.9648711 0.5799015 2.2160546 14.5996889 15.1806222
43 44 45 46 47 48
12.6861382 22.8681639 20.0698255 14.9427999 15.3060928 14.2189682
49 50 51 52 53 54
6.7520714 -2.5778829 -5.2364112 -5.3609939 -14.1118436 -12.7423671
55 56 57 58 59 60
-22.1220560 -18.7875645 -13.9563885 -20.4996889 -2.9560471 -19.0548557
> postscript(file="/var/www/html/rcomp/tmp/6vk1w1258637484.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 = 60
Frequency = 1
lag(myerror, k = 1) myerror
0 5.2718739 NA
1 6.5109633 5.2718739
2 5.9849847 6.5109633
3 -5.9824277 5.9849847
4 -4.9871852 -5.9824277
5 -0.1372460 -4.9871852
6 -3.5562749 -0.1372460
7 -4.0350229 -3.5562749
8 -9.4940971 -4.0350229
9 0.6561910 -9.4940971
10 -9.5018742 0.6561910
11 -5.3563355 -9.5018742
12 -9.7511154 -5.3563355
13 -13.6104704 -9.7511154
14 -4.5429363 -13.6104704
15 1.8636645 -4.5429363
16 -3.3872230 1.8636645
17 -4.3176329 -3.3872230
18 0.3616691 -4.3176329
19 -5.2285356 0.3616691
20 -2.1596889 -5.2285356
21 1.2216691 -2.1596889
22 -4.1396586 1.2216691
23 1.8551213 -4.1396586
24 -1.3478528 1.8551213
25 0.7125189 -1.3478528
26 3.2144612 0.7125189
27 7.2637024 3.2144612
28 7.8865630 7.2637024
29 2.0166238 7.8865630
30 12.6305235 2.0166238
31 5.1829590 12.6305235
32 5.5403490 5.1829590
33 3.6790288 5.5403490
34 1.2914871 3.6790288
35 8.3371016 1.2914871
36 -0.9249771 8.3371016
37 8.9648711 -0.9249771
38 0.5799015 8.9648711
39 2.2160546 0.5799015
40 14.5996889 2.2160546
41 15.1806222 14.5996889
42 12.6861382 15.1806222
43 22.8681639 12.6861382
44 20.0698255 22.8681639
45 14.9427999 20.0698255
46 15.3060928 14.9427999
47 14.2189682 15.3060928
48 6.7520714 14.2189682
49 -2.5778829 6.7520714
50 -5.2364112 -2.5778829
51 -5.3609939 -5.2364112
52 -14.1118436 -5.3609939
53 -12.7423671 -14.1118436
54 -22.1220560 -12.7423671
55 -18.7875645 -22.1220560
56 -13.9563885 -18.7875645
57 -20.4996889 -13.9563885
58 -2.9560471 -20.4996889
59 -19.0548557 -2.9560471
60 NA -19.0548557
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 6.5109633 5.2718739
[2,] 5.9849847 6.5109633
[3,] -5.9824277 5.9849847
[4,] -4.9871852 -5.9824277
[5,] -0.1372460 -4.9871852
[6,] -3.5562749 -0.1372460
[7,] -4.0350229 -3.5562749
[8,] -9.4940971 -4.0350229
[9,] 0.6561910 -9.4940971
[10,] -9.5018742 0.6561910
[11,] -5.3563355 -9.5018742
[12,] -9.7511154 -5.3563355
[13,] -13.6104704 -9.7511154
[14,] -4.5429363 -13.6104704
[15,] 1.8636645 -4.5429363
[16,] -3.3872230 1.8636645
[17,] -4.3176329 -3.3872230
[18,] 0.3616691 -4.3176329
[19,] -5.2285356 0.3616691
[20,] -2.1596889 -5.2285356
[21,] 1.2216691 -2.1596889
[22,] -4.1396586 1.2216691
[23,] 1.8551213 -4.1396586
[24,] -1.3478528 1.8551213
[25,] 0.7125189 -1.3478528
[26,] 3.2144612 0.7125189
[27,] 7.2637024 3.2144612
[28,] 7.8865630 7.2637024
[29,] 2.0166238 7.8865630
[30,] 12.6305235 2.0166238
[31,] 5.1829590 12.6305235
[32,] 5.5403490 5.1829590
[33,] 3.6790288 5.5403490
[34,] 1.2914871 3.6790288
[35,] 8.3371016 1.2914871
[36,] -0.9249771 8.3371016
[37,] 8.9648711 -0.9249771
[38,] 0.5799015 8.9648711
[39,] 2.2160546 0.5799015
[40,] 14.5996889 2.2160546
[41,] 15.1806222 14.5996889
[42,] 12.6861382 15.1806222
[43,] 22.8681639 12.6861382
[44,] 20.0698255 22.8681639
[45,] 14.9427999 20.0698255
[46,] 15.3060928 14.9427999
[47,] 14.2189682 15.3060928
[48,] 6.7520714 14.2189682
[49,] -2.5778829 6.7520714
[50,] -5.2364112 -2.5778829
[51,] -5.3609939 -5.2364112
[52,] -14.1118436 -5.3609939
[53,] -12.7423671 -14.1118436
[54,] -22.1220560 -12.7423671
[55,] -18.7875645 -22.1220560
[56,] -13.9563885 -18.7875645
[57,] -20.4996889 -13.9563885
[58,] -2.9560471 -20.4996889
[59,] -19.0548557 -2.9560471
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 6.5109633 5.2718739
2 5.9849847 6.5109633
3 -5.9824277 5.9849847
4 -4.9871852 -5.9824277
5 -0.1372460 -4.9871852
6 -3.5562749 -0.1372460
7 -4.0350229 -3.5562749
8 -9.4940971 -4.0350229
9 0.6561910 -9.4940971
10 -9.5018742 0.6561910
11 -5.3563355 -9.5018742
12 -9.7511154 -5.3563355
13 -13.6104704 -9.7511154
14 -4.5429363 -13.6104704
15 1.8636645 -4.5429363
16 -3.3872230 1.8636645
17 -4.3176329 -3.3872230
18 0.3616691 -4.3176329
19 -5.2285356 0.3616691
20 -2.1596889 -5.2285356
21 1.2216691 -2.1596889
22 -4.1396586 1.2216691
23 1.8551213 -4.1396586
24 -1.3478528 1.8551213
25 0.7125189 -1.3478528
26 3.2144612 0.7125189
27 7.2637024 3.2144612
28 7.8865630 7.2637024
29 2.0166238 7.8865630
30 12.6305235 2.0166238
31 5.1829590 12.6305235
32 5.5403490 5.1829590
33 3.6790288 5.5403490
34 1.2914871 3.6790288
35 8.3371016 1.2914871
36 -0.9249771 8.3371016
37 8.9648711 -0.9249771
38 0.5799015 8.9648711
39 2.2160546 0.5799015
40 14.5996889 2.2160546
41 15.1806222 14.5996889
42 12.6861382 15.1806222
43 22.8681639 12.6861382
44 20.0698255 22.8681639
45 14.9427999 20.0698255
46 15.3060928 14.9427999
47 14.2189682 15.3060928
48 6.7520714 14.2189682
49 -2.5778829 6.7520714
50 -5.2364112 -2.5778829
51 -5.3609939 -5.2364112
52 -14.1118436 -5.3609939
53 -12.7423671 -14.1118436
54 -22.1220560 -12.7423671
55 -18.7875645 -22.1220560
56 -13.9563885 -18.7875645
57 -20.4996889 -13.9563885
58 -2.9560471 -20.4996889
59 -19.0548557 -2.9560471
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/7fdt11258637484.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/8pwhg1258637484.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/9myjl1258637484.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/rcomp/tmp/10q9e11258637484.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/1140he1258637484.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/121zmb1258637484.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/13g1b11258637484.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/14af6y1258637484.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/15jzr31258637484.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/16y8p41258637484.tab")
+ }
>
> system("convert tmp/19zj41258637484.ps tmp/19zj41258637484.png")
> system("convert tmp/2gkm41258637484.ps tmp/2gkm41258637484.png")
> system("convert tmp/3ag0g1258637484.ps tmp/3ag0g1258637484.png")
> system("convert tmp/4cp4j1258637484.ps tmp/4cp4j1258637484.png")
> system("convert tmp/50xb01258637484.ps tmp/50xb01258637484.png")
> system("convert tmp/6vk1w1258637484.ps tmp/6vk1w1258637484.png")
> system("convert tmp/7fdt11258637484.ps tmp/7fdt11258637484.png")
> system("convert tmp/8pwhg1258637484.ps tmp/8pwhg1258637484.png")
> system("convert tmp/9myjl1258637484.ps tmp/9myjl1258637484.png")
> system("convert tmp/10q9e11258637484.ps tmp/10q9e11258637484.png")
>
>
> proc.time()
user system elapsed
2.215 1.651 2.880