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(108.2,108.5,108.8,112.3,110.2,116.6,109.5,115.5,109.5,120.1,116,132.9,111.2,128.1,112.1,129.3,114,132.5,119.1,131,114.1,124.9,115.1,120.8,115.4,122,110.8,122.1,116,127.4,119.2,135.2,126.5,137.3,127.8,135,131.3,136,140.3,138.4,137.3,134.7,143,138.4,134.5,133.9,139.9,133.6,159.3,141.2,170.4,151.8,175,155.4,175.8,156.6,180.9,161.6,180.3,160.7,169.6,156,172.3,159.5,184.8,168.7,177.7,169.9,184.6,169.9,211.4,185.9,215.3,190.8,215.9,195.8,244.7,211.9,259.3,227.1,289,251.3,310.9,256.7,321,251.9,315.1,251.2,333.2,270.3,314.1,267.2,284.7,243,273.9,229.9,216,187.2,196.4,178.2,190.9,175.2,206.4,192.4,196.3,187,199.5,184,198.9,194.1,214.4,212.7,214.2,217.5,187.6,200.5,180.6,205.9,172.2,196.5),dim=c(2,60),dimnames=list(c('Y','X'),1:60))
> y <- array(NA,dim=c(2,60),dimnames=list(c('Y','X'),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
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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 108.2 108.5 1 0 0 0 0 0 0 0 0 0 0 1
2 108.8 112.3 0 1 0 0 0 0 0 0 0 0 0 2
3 110.2 116.6 0 0 1 0 0 0 0 0 0 0 0 3
4 109.5 115.5 0 0 0 1 0 0 0 0 0 0 0 4
5 109.5 120.1 0 0 0 0 1 0 0 0 0 0 0 5
6 116.0 132.9 0 0 0 0 0 1 0 0 0 0 0 6
7 111.2 128.1 0 0 0 0 0 0 1 0 0 0 0 7
8 112.1 129.3 0 0 0 0 0 0 0 1 0 0 0 8
9 114.0 132.5 0 0 0 0 0 0 0 0 1 0 0 9
10 119.1 131.0 0 0 0 0 0 0 0 0 0 1 0 10
11 114.1 124.9 0 0 0 0 0 0 0 0 0 0 1 11
12 115.1 120.8 0 0 0 0 0 0 0 0 0 0 0 12
13 115.4 122.0 1 0 0 0 0 0 0 0 0 0 0 13
14 110.8 122.1 0 1 0 0 0 0 0 0 0 0 0 14
15 116.0 127.4 0 0 1 0 0 0 0 0 0 0 0 15
16 119.2 135.2 0 0 0 1 0 0 0 0 0 0 0 16
17 126.5 137.3 0 0 0 0 1 0 0 0 0 0 0 17
18 127.8 135.0 0 0 0 0 0 1 0 0 0 0 0 18
19 131.3 136.0 0 0 0 0 0 0 1 0 0 0 0 19
20 140.3 138.4 0 0 0 0 0 0 0 1 0 0 0 20
21 137.3 134.7 0 0 0 0 0 0 0 0 1 0 0 21
22 143.0 138.4 0 0 0 0 0 0 0 0 0 1 0 22
23 134.5 133.9 0 0 0 0 0 0 0 0 0 0 1 23
24 139.9 133.6 0 0 0 0 0 0 0 0 0 0 0 24
25 159.3 141.2 1 0 0 0 0 0 0 0 0 0 0 25
26 170.4 151.8 0 1 0 0 0 0 0 0 0 0 0 26
27 175.0 155.4 0 0 1 0 0 0 0 0 0 0 0 27
28 175.8 156.6 0 0 0 1 0 0 0 0 0 0 0 28
29 180.9 161.6 0 0 0 0 1 0 0 0 0 0 0 29
30 180.3 160.7 0 0 0 0 0 1 0 0 0 0 0 30
31 169.6 156.0 0 0 0 0 0 0 1 0 0 0 0 31
32 172.3 159.5 0 0 0 0 0 0 0 1 0 0 0 32
33 184.8 168.7 0 0 0 0 0 0 0 0 1 0 0 33
34 177.7 169.9 0 0 0 0 0 0 0 0 0 1 0 34
35 184.6 169.9 0 0 0 0 0 0 0 0 0 0 1 35
36 211.4 185.9 0 0 0 0 0 0 0 0 0 0 0 36
37 215.3 190.8 1 0 0 0 0 0 0 0 0 0 0 37
38 215.9 195.8 0 1 0 0 0 0 0 0 0 0 0 38
39 244.7 211.9 0 0 1 0 0 0 0 0 0 0 0 39
40 259.3 227.1 0 0 0 1 0 0 0 0 0 0 0 40
41 289.0 251.3 0 0 0 0 1 0 0 0 0 0 0 41
42 310.9 256.7 0 0 0 0 0 1 0 0 0 0 0 42
43 321.0 251.9 0 0 0 0 0 0 1 0 0 0 0 43
44 315.1 251.2 0 0 0 0 0 0 0 1 0 0 0 44
45 333.2 270.3 0 0 0 0 0 0 0 0 1 0 0 45
46 314.1 267.2 0 0 0 0 0 0 0 0 0 1 0 46
47 284.7 243.0 0 0 0 0 0 0 0 0 0 0 1 47
48 273.9 229.9 0 0 0 0 0 0 0 0 0 0 0 48
49 216.0 187.2 1 0 0 0 0 0 0 0 0 0 0 49
50 196.4 178.2 0 1 0 0 0 0 0 0 0 0 0 50
51 190.9 175.2 0 0 1 0 0 0 0 0 0 0 0 51
52 206.4 192.4 0 0 0 1 0 0 0 0 0 0 0 52
53 196.3 187.0 0 0 0 0 1 0 0 0 0 0 0 53
54 199.5 184.0 0 0 0 0 0 1 0 0 0 0 0 54
55 198.9 194.1 0 0 0 0 0 0 1 0 0 0 0 55
56 214.4 212.7 0 0 0 0 0 0 0 1 0 0 0 56
57 214.2 217.5 0 0 0 0 0 0 0 0 1 0 0 57
58 187.6 200.5 0 0 0 0 0 0 0 0 0 1 0 58
59 180.6 205.9 0 0 0 0 0 0 0 0 0 0 1 59
60 172.2 196.5 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) X M1 M2 M3 M4
-73.15102 1.59969 11.16068 6.02243 5.10915 -0.50327
M5 M6 M7 M8 M9 M10
-3.26029 -0.03845 1.08646 -1.87090 -5.83980 -8.29572
M11 t
-6.88843 -0.60110
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-33.2722 -6.7181 0.8721 10.5008 16.2736
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -73.15102 10.22851 -7.152 5.4e-09 ***
X 1.59969 0.07222 22.150 < 2e-16 ***
M1 11.16068 8.89208 1.255 0.21577
M2 6.02243 8.87878 0.678 0.50098
M3 5.10915 8.86848 0.576 0.56735
M4 -0.50327 8.87719 -0.057 0.95504
M5 -3.26029 8.89268 -0.367 0.71558
M6 -0.03845 8.88698 -0.004 0.99657
M7 1.08646 8.86167 0.123 0.90296
M8 -1.87090 8.87611 -0.211 0.83399
M9 -5.83980 8.91184 -0.655 0.51555
M10 -8.29572 8.86147 -0.936 0.35408
M11 -6.88843 8.82107 -0.781 0.43886
t -0.60110 0.18286 -3.287 0.00194 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.94 on 46 degrees of freedom
Multiple R-squared: 0.9623, Adjusted R-squared: 0.9517
F-statistic: 90.41 on 13 and 46 DF, p-value: < 2.2e-16
> 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.0182214542 0.036442908 0.981778546
[2,] 0.0223383755 0.044676751 0.977661625
[3,] 0.0340856820 0.068171364 0.965914318
[4,] 0.0832368997 0.166473799 0.916763100
[5,] 0.0521837516 0.104367503 0.947816248
[6,] 0.0340875532 0.068175106 0.965912447
[7,] 0.0175753555 0.035150711 0.982424645
[8,] 0.0136325893 0.027265179 0.986367411
[9,] 0.0541206130 0.108241226 0.945879387
[10,] 0.0790382884 0.158076577 0.920961712
[11,] 0.0698140093 0.139628019 0.930185991
[12,] 0.0538865463 0.107773093 0.946113454
[13,] 0.0361243051 0.072248610 0.963875695
[14,] 0.0275220949 0.055044190 0.972477905
[15,] 0.0176591055 0.035318211 0.982340895
[16,] 0.0096869954 0.019373991 0.990313005
[17,] 0.0047269325 0.009453865 0.995273068
[18,] 0.0036427205 0.007285441 0.996357279
[19,] 0.0018756181 0.003751236 0.998124382
[20,] 0.0009370807 0.001874161 0.999062919
[21,] 0.0009921792 0.001984358 0.999007821
[22,] 0.0033549916 0.006709983 0.996645008
[23,] 0.0117907893 0.023581579 0.988209211
[24,] 0.3353477740 0.670695548 0.664652226
[25,] 0.7642603165 0.471479367 0.235739684
[26,] 0.9936064854 0.012787029 0.006393515
[27,] 0.9827755588 0.034448882 0.017224441
> postscript(file="/var/www/html/rcomp/tmp/1jx171258729051.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/2az7e1258729051.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/3fbv41258729051.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/45thv1258729051.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/5hanv1258729051.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
-2.7751843 -2.5146613 -6.4789566 0.7942251 -3.2062364 -19.8030366
7 8 9 10 11 12
-17.4483165 -14.9094857 -13.5585071 -3.0019410 0.9499913 2.2214006
13 14 15 16 17 18
-9.9578135 -8.9784288 -10.7424165 -13.8064965 -6.5077272 -4.1491732
19 20 21 22 23 24
-2.7726687 5.9465314 13.4353870 16.2735530 14.1659776 13.7585560
25 26 27 28 29 30
10.4413110 10.3239262 10.6794154 15.7733048 16.2329664 14.4519511
31 32 33 34 35 36
10.7467020 11.4062405 13.7590650 7.7964618 13.8902709 8.8078643
37 38 39 40 41 42
-5.6902114 -7.3493192 -2.7899841 -6.2917873 -11.9462185 -1.3052955
43 44 45 46 47 48
15.9494246 14.7276709 6.8435413 -4.2403848 4.2659787 8.1346190
49 50 51 52 53 54
7.9818983 8.5184831 9.3319418 3.5307539 5.4272157 10.8055543
55 56 57 58 59 60
-6.4751413 -17.1709571 -20.4794863 -16.8276890 -33.2722185 -32.9224399
> postscript(file="/var/www/html/rcomp/tmp/6vy731258729051.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 -2.7751843 NA
1 -2.5146613 -2.7751843
2 -6.4789566 -2.5146613
3 0.7942251 -6.4789566
4 -3.2062364 0.7942251
5 -19.8030366 -3.2062364
6 -17.4483165 -19.8030366
7 -14.9094857 -17.4483165
8 -13.5585071 -14.9094857
9 -3.0019410 -13.5585071
10 0.9499913 -3.0019410
11 2.2214006 0.9499913
12 -9.9578135 2.2214006
13 -8.9784288 -9.9578135
14 -10.7424165 -8.9784288
15 -13.8064965 -10.7424165
16 -6.5077272 -13.8064965
17 -4.1491732 -6.5077272
18 -2.7726687 -4.1491732
19 5.9465314 -2.7726687
20 13.4353870 5.9465314
21 16.2735530 13.4353870
22 14.1659776 16.2735530
23 13.7585560 14.1659776
24 10.4413110 13.7585560
25 10.3239262 10.4413110
26 10.6794154 10.3239262
27 15.7733048 10.6794154
28 16.2329664 15.7733048
29 14.4519511 16.2329664
30 10.7467020 14.4519511
31 11.4062405 10.7467020
32 13.7590650 11.4062405
33 7.7964618 13.7590650
34 13.8902709 7.7964618
35 8.8078643 13.8902709
36 -5.6902114 8.8078643
37 -7.3493192 -5.6902114
38 -2.7899841 -7.3493192
39 -6.2917873 -2.7899841
40 -11.9462185 -6.2917873
41 -1.3052955 -11.9462185
42 15.9494246 -1.3052955
43 14.7276709 15.9494246
44 6.8435413 14.7276709
45 -4.2403848 6.8435413
46 4.2659787 -4.2403848
47 8.1346190 4.2659787
48 7.9818983 8.1346190
49 8.5184831 7.9818983
50 9.3319418 8.5184831
51 3.5307539 9.3319418
52 5.4272157 3.5307539
53 10.8055543 5.4272157
54 -6.4751413 10.8055543
55 -17.1709571 -6.4751413
56 -20.4794863 -17.1709571
57 -16.8276890 -20.4794863
58 -33.2722185 -16.8276890
59 -32.9224399 -33.2722185
60 NA -32.9224399
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2.5146613 -2.7751843
[2,] -6.4789566 -2.5146613
[3,] 0.7942251 -6.4789566
[4,] -3.2062364 0.7942251
[5,] -19.8030366 -3.2062364
[6,] -17.4483165 -19.8030366
[7,] -14.9094857 -17.4483165
[8,] -13.5585071 -14.9094857
[9,] -3.0019410 -13.5585071
[10,] 0.9499913 -3.0019410
[11,] 2.2214006 0.9499913
[12,] -9.9578135 2.2214006
[13,] -8.9784288 -9.9578135
[14,] -10.7424165 -8.9784288
[15,] -13.8064965 -10.7424165
[16,] -6.5077272 -13.8064965
[17,] -4.1491732 -6.5077272
[18,] -2.7726687 -4.1491732
[19,] 5.9465314 -2.7726687
[20,] 13.4353870 5.9465314
[21,] 16.2735530 13.4353870
[22,] 14.1659776 16.2735530
[23,] 13.7585560 14.1659776
[24,] 10.4413110 13.7585560
[25,] 10.3239262 10.4413110
[26,] 10.6794154 10.3239262
[27,] 15.7733048 10.6794154
[28,] 16.2329664 15.7733048
[29,] 14.4519511 16.2329664
[30,] 10.7467020 14.4519511
[31,] 11.4062405 10.7467020
[32,] 13.7590650 11.4062405
[33,] 7.7964618 13.7590650
[34,] 13.8902709 7.7964618
[35,] 8.8078643 13.8902709
[36,] -5.6902114 8.8078643
[37,] -7.3493192 -5.6902114
[38,] -2.7899841 -7.3493192
[39,] -6.2917873 -2.7899841
[40,] -11.9462185 -6.2917873
[41,] -1.3052955 -11.9462185
[42,] 15.9494246 -1.3052955
[43,] 14.7276709 15.9494246
[44,] 6.8435413 14.7276709
[45,] -4.2403848 6.8435413
[46,] 4.2659787 -4.2403848
[47,] 8.1346190 4.2659787
[48,] 7.9818983 8.1346190
[49,] 8.5184831 7.9818983
[50,] 9.3319418 8.5184831
[51,] 3.5307539 9.3319418
[52,] 5.4272157 3.5307539
[53,] 10.8055543 5.4272157
[54,] -6.4751413 10.8055543
[55,] -17.1709571 -6.4751413
[56,] -20.4794863 -17.1709571
[57,] -16.8276890 -20.4794863
[58,] -33.2722185 -16.8276890
[59,] -32.9224399 -33.2722185
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2.5146613 -2.7751843
2 -6.4789566 -2.5146613
3 0.7942251 -6.4789566
4 -3.2062364 0.7942251
5 -19.8030366 -3.2062364
6 -17.4483165 -19.8030366
7 -14.9094857 -17.4483165
8 -13.5585071 -14.9094857
9 -3.0019410 -13.5585071
10 0.9499913 -3.0019410
11 2.2214006 0.9499913
12 -9.9578135 2.2214006
13 -8.9784288 -9.9578135
14 -10.7424165 -8.9784288
15 -13.8064965 -10.7424165
16 -6.5077272 -13.8064965
17 -4.1491732 -6.5077272
18 -2.7726687 -4.1491732
19 5.9465314 -2.7726687
20 13.4353870 5.9465314
21 16.2735530 13.4353870
22 14.1659776 16.2735530
23 13.7585560 14.1659776
24 10.4413110 13.7585560
25 10.3239262 10.4413110
26 10.6794154 10.3239262
27 15.7733048 10.6794154
28 16.2329664 15.7733048
29 14.4519511 16.2329664
30 10.7467020 14.4519511
31 11.4062405 10.7467020
32 13.7590650 11.4062405
33 7.7964618 13.7590650
34 13.8902709 7.7964618
35 8.8078643 13.8902709
36 -5.6902114 8.8078643
37 -7.3493192 -5.6902114
38 -2.7899841 -7.3493192
39 -6.2917873 -2.7899841
40 -11.9462185 -6.2917873
41 -1.3052955 -11.9462185
42 15.9494246 -1.3052955
43 14.7276709 15.9494246
44 6.8435413 14.7276709
45 -4.2403848 6.8435413
46 4.2659787 -4.2403848
47 8.1346190 4.2659787
48 7.9818983 8.1346190
49 8.5184831 7.9818983
50 9.3319418 8.5184831
51 3.5307539 9.3319418
52 5.4272157 3.5307539
53 10.8055543 5.4272157
54 -6.4751413 10.8055543
55 -17.1709571 -6.4751413
56 -20.4794863 -17.1709571
57 -16.8276890 -20.4794863
58 -33.2722185 -16.8276890
59 -32.9224399 -33.2722185
> 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/76neb1258729051.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/8lkx31258729051.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/9z4c71258729051.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/105qf01258729051.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/11qimw1258729051.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/12m91u1258729051.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/13uzim1258729051.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/14wyir1258729051.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/15gp7b1258729051.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/166z8k1258729051.tab")
+ }
>
> system("convert tmp/1jx171258729051.ps tmp/1jx171258729051.png")
> system("convert tmp/2az7e1258729051.ps tmp/2az7e1258729051.png")
> system("convert tmp/3fbv41258729051.ps tmp/3fbv41258729051.png")
> system("convert tmp/45thv1258729051.ps tmp/45thv1258729051.png")
> system("convert tmp/5hanv1258729051.ps tmp/5hanv1258729051.png")
> system("convert tmp/6vy731258729051.ps tmp/6vy731258729051.png")
> system("convert tmp/76neb1258729051.ps tmp/76neb1258729051.png")
> system("convert tmp/8lkx31258729051.ps tmp/8lkx31258729051.png")
> system("convert tmp/9z4c71258729051.ps tmp/9z4c71258729051.png")
> system("convert tmp/105qf01258729051.ps tmp/105qf01258729051.png")
>
>
> proc.time()
user system elapsed
2.392 1.546 2.784