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 = '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
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 108.2 108.5 1 0 0 0 0 0 0 0 0 0 0
2 108.8 112.3 0 1 0 0 0 0 0 0 0 0 0
3 110.2 116.6 0 0 1 0 0 0 0 0 0 0 0
4 109.5 115.5 0 0 0 1 0 0 0 0 0 0 0
5 109.5 120.1 0 0 0 0 1 0 0 0 0 0 0
6 116.0 132.9 0 0 0 0 0 1 0 0 0 0 0
7 111.2 128.1 0 0 0 0 0 0 1 0 0 0 0
8 112.1 129.3 0 0 0 0 0 0 0 1 0 0 0
9 114.0 132.5 0 0 0 0 0 0 0 0 1 0 0
10 119.1 131.0 0 0 0 0 0 0 0 0 0 1 0
11 114.1 124.9 0 0 0 0 0 0 0 0 0 0 1
12 115.1 120.8 0 0 0 0 0 0 0 0 0 0 0
13 115.4 122.0 1 0 0 0 0 0 0 0 0 0 0
14 110.8 122.1 0 1 0 0 0 0 0 0 0 0 0
15 116.0 127.4 0 0 1 0 0 0 0 0 0 0 0
16 119.2 135.2 0 0 0 1 0 0 0 0 0 0 0
17 126.5 137.3 0 0 0 0 1 0 0 0 0 0 0
18 127.8 135.0 0 0 0 0 0 1 0 0 0 0 0
19 131.3 136.0 0 0 0 0 0 0 1 0 0 0 0
20 140.3 138.4 0 0 0 0 0 0 0 1 0 0 0
21 137.3 134.7 0 0 0 0 0 0 0 0 1 0 0
22 143.0 138.4 0 0 0 0 0 0 0 0 0 1 0
23 134.5 133.9 0 0 0 0 0 0 0 0 0 0 1
24 139.9 133.6 0 0 0 0 0 0 0 0 0 0 0
25 159.3 141.2 1 0 0 0 0 0 0 0 0 0 0
26 170.4 151.8 0 1 0 0 0 0 0 0 0 0 0
27 175.0 155.4 0 0 1 0 0 0 0 0 0 0 0
28 175.8 156.6 0 0 0 1 0 0 0 0 0 0 0
29 180.9 161.6 0 0 0 0 1 0 0 0 0 0 0
30 180.3 160.7 0 0 0 0 0 1 0 0 0 0 0
31 169.6 156.0 0 0 0 0 0 0 1 0 0 0 0
32 172.3 159.5 0 0 0 0 0 0 0 1 0 0 0
33 184.8 168.7 0 0 0 0 0 0 0 0 1 0 0
34 177.7 169.9 0 0 0 0 0 0 0 0 0 1 0
35 184.6 169.9 0 0 0 0 0 0 0 0 0 0 1
36 211.4 185.9 0 0 0 0 0 0 0 0 0 0 0
37 215.3 190.8 1 0 0 0 0 0 0 0 0 0 0
38 215.9 195.8 0 1 0 0 0 0 0 0 0 0 0
39 244.7 211.9 0 0 1 0 0 0 0 0 0 0 0
40 259.3 227.1 0 0 0 1 0 0 0 0 0 0 0
41 289.0 251.3 0 0 0 0 1 0 0 0 0 0 0
42 310.9 256.7 0 0 0 0 0 1 0 0 0 0 0
43 321.0 251.9 0 0 0 0 0 0 1 0 0 0 0
44 315.1 251.2 0 0 0 0 0 0 0 1 0 0 0
45 333.2 270.3 0 0 0 0 0 0 0 0 1 0 0
46 314.1 267.2 0 0 0 0 0 0 0 0 0 1 0
47 284.7 243.0 0 0 0 0 0 0 0 0 0 0 1
48 273.9 229.9 0 0 0 0 0 0 0 0 0 0 0
49 216.0 187.2 1 0 0 0 0 0 0 0 0 0 0
50 196.4 178.2 0 1 0 0 0 0 0 0 0 0 0
51 190.9 175.2 0 0 1 0 0 0 0 0 0 0 0
52 206.4 192.4 0 0 0 1 0 0 0 0 0 0 0
53 196.3 187.0 0 0 0 0 1 0 0 0 0 0 0
54 199.5 184.0 0 0 0 0 0 1 0 0 0 0 0
55 198.9 194.1 0 0 0 0 0 0 1 0 0 0 0
56 214.4 212.7 0 0 0 0 0 0 0 1 0 0 0
57 214.2 217.5 0 0 0 0 0 0 0 0 1 0 0
58 187.6 200.5 0 0 0 0 0 0 0 0 0 1 0
59 180.6 205.9 0 0 0 0 0 0 0 0 0 0 1
60 172.2 196.5 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) X M1 M2 M3 M4
-61.2630 1.4063 13.2467 7.9136 7.4166 2.7620
M5 M6 M7 M8 M9 M10
0.5838 3.6687 4.0688 1.4774 -1.8315 -5.5345
M11
-5.8657
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-42.869 -6.020 1.481 10.286 23.955
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -61.26304 10.51875 -5.824 4.96e-07 ***
X 1.40627 0.04604 30.546 < 2e-16 ***
M1 13.24674 9.75086 1.359 0.181
M2 7.91358 9.74066 0.812 0.421
M3 7.41659 9.71926 0.763 0.449
M4 2.76204 9.69813 0.285 0.777
M5 0.58379 9.69156 0.060 0.952
M6 3.66874 9.69120 0.379 0.707
M7 4.06875 9.69117 0.420 0.677
M8 1.47740 9.69377 0.152 0.880
M9 -1.83149 9.70537 -0.189 0.851
M10 -5.53455 9.69827 -0.571 0.571
M11 -5.86567 9.69169 -0.605 0.548
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15.32 on 47 degrees of freedom
Multiple R-squared: 0.9535, Adjusted R-squared: 0.9416
F-statistic: 80.29 on 12 and 47 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,] 4.073580e-04 0.0008147160 0.9995926
[2,] 1.318369e-03 0.0026367389 0.9986816
[3,] 1.987210e-03 0.0039744208 0.9980128
[4,] 3.852512e-03 0.0077050238 0.9961475
[5,] 1.041306e-02 0.0208261283 0.9895869
[6,] 1.601008e-02 0.0320201562 0.9839899
[7,] 1.426647e-02 0.0285329343 0.9857335
[8,] 8.845690e-03 0.0176913798 0.9911543
[9,] 5.447875e-03 0.0108957504 0.9945521
[10,] 4.098424e-03 0.0081968474 0.9959016
[11,] 2.332965e-03 0.0046659305 0.9976670
[12,] 1.262182e-03 0.0025243647 0.9987378
[13,] 7.932850e-04 0.0015865701 0.9992067
[14,] 4.605104e-04 0.0009210208 0.9995395
[15,] 3.034989e-04 0.0006069979 0.9996965
[16,] 1.342535e-04 0.0002685070 0.9998657
[17,] 6.087513e-05 0.0001217503 0.9999391
[18,] 4.559755e-05 0.0000911951 0.9999544
[19,] 9.932035e-05 0.0001986407 0.9999007
[20,] 5.566841e-04 0.0011133683 0.9994433
[21,] 2.829434e-03 0.0056588688 0.9971706
[22,] 3.554565e-03 0.0071091297 0.9964454
[23,] 3.737798e-03 0.0074755969 0.9962622
[24,] 2.482009e-03 0.0049640184 0.9975180
[25,] 1.650447e-03 0.0033008940 0.9983496
[26,] 4.907339e-03 0.0098146785 0.9950927
[27,] 6.106370e-02 0.1221274066 0.9389363
[28,] 5.552313e-02 0.1110462664 0.9444769
[29,] 3.913799e-02 0.0782759773 0.9608620
> postscript(file="/var/www/html/rcomp/tmp/1tntb1258727779.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/2oigp1258727779.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/3854s1258727779.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/42wda1258727779.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/5q45i1258727779.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
3.6358756 4.2252147 0.0752350 5.5766785 1.2860852 -13.2991345
7 8 9 10 11 12
-11.7490466 -9.9452163 -9.2363962 1.6760650 5.5854446 6.4854852
13 14 15 16 17 18
-8.1487846 -7.5562424 -9.3124932 -12.4268627 -5.9017782 -4.4523039
19 20 21 22 23 24
-2.7585885 5.4577164 10.9698073 15.1696586 13.3290044 13.2852147
25 26 27 28 29 30
8.7508097 10.2775051 10.3119151 14.0789351 14.3258333 11.9065281
31 32 33 34 35 36
7.4159889 7.7853956 10.6565889 5.5721180 12.8032438 11.2372346
37 38 39 40 41 42
-5.0002384 -6.0984247 0.5575963 -1.5631796 -3.7166870 7.5044996
43 44 45 46 47 48
23.9545875 21.6303329 16.1794421 5.1419370 10.1048241 11.8613049
49 50 51 52 53 54
0.7623377 -0.8480528 -1.6322532 -5.6655714 -5.9934534 -1.6595893
55 56 57 58 59 60
-16.8629412 -24.9282286 -28.5694422 -27.5597786 -41.8225169 -42.8692394
> postscript(file="/var/www/html/rcomp/tmp/6qgi41258727779.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 3.6358756 NA
1 4.2252147 3.6358756
2 0.0752350 4.2252147
3 5.5766785 0.0752350
4 1.2860852 5.5766785
5 -13.2991345 1.2860852
6 -11.7490466 -13.2991345
7 -9.9452163 -11.7490466
8 -9.2363962 -9.9452163
9 1.6760650 -9.2363962
10 5.5854446 1.6760650
11 6.4854852 5.5854446
12 -8.1487846 6.4854852
13 -7.5562424 -8.1487846
14 -9.3124932 -7.5562424
15 -12.4268627 -9.3124932
16 -5.9017782 -12.4268627
17 -4.4523039 -5.9017782
18 -2.7585885 -4.4523039
19 5.4577164 -2.7585885
20 10.9698073 5.4577164
21 15.1696586 10.9698073
22 13.3290044 15.1696586
23 13.2852147 13.3290044
24 8.7508097 13.2852147
25 10.2775051 8.7508097
26 10.3119151 10.2775051
27 14.0789351 10.3119151
28 14.3258333 14.0789351
29 11.9065281 14.3258333
30 7.4159889 11.9065281
31 7.7853956 7.4159889
32 10.6565889 7.7853956
33 5.5721180 10.6565889
34 12.8032438 5.5721180
35 11.2372346 12.8032438
36 -5.0002384 11.2372346
37 -6.0984247 -5.0002384
38 0.5575963 -6.0984247
39 -1.5631796 0.5575963
40 -3.7166870 -1.5631796
41 7.5044996 -3.7166870
42 23.9545875 7.5044996
43 21.6303329 23.9545875
44 16.1794421 21.6303329
45 5.1419370 16.1794421
46 10.1048241 5.1419370
47 11.8613049 10.1048241
48 0.7623377 11.8613049
49 -0.8480528 0.7623377
50 -1.6322532 -0.8480528
51 -5.6655714 -1.6322532
52 -5.9934534 -5.6655714
53 -1.6595893 -5.9934534
54 -16.8629412 -1.6595893
55 -24.9282286 -16.8629412
56 -28.5694422 -24.9282286
57 -27.5597786 -28.5694422
58 -41.8225169 -27.5597786
59 -42.8692394 -41.8225169
60 NA -42.8692394
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 4.2252147 3.6358756
[2,] 0.0752350 4.2252147
[3,] 5.5766785 0.0752350
[4,] 1.2860852 5.5766785
[5,] -13.2991345 1.2860852
[6,] -11.7490466 -13.2991345
[7,] -9.9452163 -11.7490466
[8,] -9.2363962 -9.9452163
[9,] 1.6760650 -9.2363962
[10,] 5.5854446 1.6760650
[11,] 6.4854852 5.5854446
[12,] -8.1487846 6.4854852
[13,] -7.5562424 -8.1487846
[14,] -9.3124932 -7.5562424
[15,] -12.4268627 -9.3124932
[16,] -5.9017782 -12.4268627
[17,] -4.4523039 -5.9017782
[18,] -2.7585885 -4.4523039
[19,] 5.4577164 -2.7585885
[20,] 10.9698073 5.4577164
[21,] 15.1696586 10.9698073
[22,] 13.3290044 15.1696586
[23,] 13.2852147 13.3290044
[24,] 8.7508097 13.2852147
[25,] 10.2775051 8.7508097
[26,] 10.3119151 10.2775051
[27,] 14.0789351 10.3119151
[28,] 14.3258333 14.0789351
[29,] 11.9065281 14.3258333
[30,] 7.4159889 11.9065281
[31,] 7.7853956 7.4159889
[32,] 10.6565889 7.7853956
[33,] 5.5721180 10.6565889
[34,] 12.8032438 5.5721180
[35,] 11.2372346 12.8032438
[36,] -5.0002384 11.2372346
[37,] -6.0984247 -5.0002384
[38,] 0.5575963 -6.0984247
[39,] -1.5631796 0.5575963
[40,] -3.7166870 -1.5631796
[41,] 7.5044996 -3.7166870
[42,] 23.9545875 7.5044996
[43,] 21.6303329 23.9545875
[44,] 16.1794421 21.6303329
[45,] 5.1419370 16.1794421
[46,] 10.1048241 5.1419370
[47,] 11.8613049 10.1048241
[48,] 0.7623377 11.8613049
[49,] -0.8480528 0.7623377
[50,] -1.6322532 -0.8480528
[51,] -5.6655714 -1.6322532
[52,] -5.9934534 -5.6655714
[53,] -1.6595893 -5.9934534
[54,] -16.8629412 -1.6595893
[55,] -24.9282286 -16.8629412
[56,] -28.5694422 -24.9282286
[57,] -27.5597786 -28.5694422
[58,] -41.8225169 -27.5597786
[59,] -42.8692394 -41.8225169
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 4.2252147 3.6358756
2 0.0752350 4.2252147
3 5.5766785 0.0752350
4 1.2860852 5.5766785
5 -13.2991345 1.2860852
6 -11.7490466 -13.2991345
7 -9.9452163 -11.7490466
8 -9.2363962 -9.9452163
9 1.6760650 -9.2363962
10 5.5854446 1.6760650
11 6.4854852 5.5854446
12 -8.1487846 6.4854852
13 -7.5562424 -8.1487846
14 -9.3124932 -7.5562424
15 -12.4268627 -9.3124932
16 -5.9017782 -12.4268627
17 -4.4523039 -5.9017782
18 -2.7585885 -4.4523039
19 5.4577164 -2.7585885
20 10.9698073 5.4577164
21 15.1696586 10.9698073
22 13.3290044 15.1696586
23 13.2852147 13.3290044
24 8.7508097 13.2852147
25 10.2775051 8.7508097
26 10.3119151 10.2775051
27 14.0789351 10.3119151
28 14.3258333 14.0789351
29 11.9065281 14.3258333
30 7.4159889 11.9065281
31 7.7853956 7.4159889
32 10.6565889 7.7853956
33 5.5721180 10.6565889
34 12.8032438 5.5721180
35 11.2372346 12.8032438
36 -5.0002384 11.2372346
37 -6.0984247 -5.0002384
38 0.5575963 -6.0984247
39 -1.5631796 0.5575963
40 -3.7166870 -1.5631796
41 7.5044996 -3.7166870
42 23.9545875 7.5044996
43 21.6303329 23.9545875
44 16.1794421 21.6303329
45 5.1419370 16.1794421
46 10.1048241 5.1419370
47 11.8613049 10.1048241
48 0.7623377 11.8613049
49 -0.8480528 0.7623377
50 -1.6322532 -0.8480528
51 -5.6655714 -1.6322532
52 -5.9934534 -5.6655714
53 -1.6595893 -5.9934534
54 -16.8629412 -1.6595893
55 -24.9282286 -16.8629412
56 -28.5694422 -24.9282286
57 -27.5597786 -28.5694422
58 -41.8225169 -27.5597786
59 -42.8692394 -41.8225169
> 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/7w5eo1258727779.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/8fpp81258727779.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/9b79x1258727779.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/10p2wv1258727779.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/115f9k1258727779.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/12gzmk1258727779.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/13kbbn1258727779.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/14gm931258727779.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/154f411258727779.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/16gl9s1258727779.tab")
+ }
>
> system("convert tmp/1tntb1258727779.ps tmp/1tntb1258727779.png")
> system("convert tmp/2oigp1258727779.ps tmp/2oigp1258727779.png")
> system("convert tmp/3854s1258727779.ps tmp/3854s1258727779.png")
> system("convert tmp/42wda1258727779.ps tmp/42wda1258727779.png")
> system("convert tmp/5q45i1258727779.ps tmp/5q45i1258727779.png")
> system("convert tmp/6qgi41258727779.ps tmp/6qgi41258727779.png")
> system("convert tmp/7w5eo1258727779.ps tmp/7w5eo1258727779.png")
> system("convert tmp/8fpp81258727779.ps tmp/8fpp81258727779.png")
> system("convert tmp/9b79x1258727779.ps tmp/9b79x1258727779.png")
> system("convert tmp/10p2wv1258727779.ps tmp/10p2wv1258727779.png")
>
>
> proc.time()
user system elapsed
2.478 1.601 5.021