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(695,0,638,0,762,0,635,0,721,0,854,0,418,0,367,0,824,0,687,0,601,0,676,0,740,0,691,0,683,0,594,0,729,0,731,0,386,0,331,0,707,0,715,0,657,0,653,0,642,0,643,0,718,0,654,0,632,0,731,0,392,1,344,1,792,1,852,1,649,1,629,1,685,1,617,1,715,1,715,1,629,1,916,1,531,1,357,1,917,1,828,1,708,1,858,1,775,1,785,1,1006,1,789,1,734,1,906,1,532,1,387,1,991,1,841,1,892,1,782,1,813,1,793,1,978,1,775,1,797,1,946,1,594,1,438,1,1022,1,868,1),dim=c(2,70),dimnames=list(c('Y','X'),1:70))
> y <- array(NA,dim=c(2,70),dimnames=list(c('Y','X'),1:70))
> 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 695 0 1 0 0 0 0 0 0 0 0 0 0
2 638 0 0 1 0 0 0 0 0 0 0 0 0
3 762 0 0 0 1 0 0 0 0 0 0 0 0
4 635 0 0 0 0 1 0 0 0 0 0 0 0
5 721 0 0 0 0 0 1 0 0 0 0 0 0
6 854 0 0 0 0 0 0 1 0 0 0 0 0
7 418 0 0 0 0 0 0 0 1 0 0 0 0
8 367 0 0 0 0 0 0 0 0 1 0 0 0
9 824 0 0 0 0 0 0 0 0 0 1 0 0
10 687 0 0 0 0 0 0 0 0 0 0 1 0
11 601 0 0 0 0 0 0 0 0 0 0 0 1
12 676 0 0 0 0 0 0 0 0 0 0 0 0
13 740 0 1 0 0 0 0 0 0 0 0 0 0
14 691 0 0 1 0 0 0 0 0 0 0 0 0
15 683 0 0 0 1 0 0 0 0 0 0 0 0
16 594 0 0 0 0 1 0 0 0 0 0 0 0
17 729 0 0 0 0 0 1 0 0 0 0 0 0
18 731 0 0 0 0 0 0 1 0 0 0 0 0
19 386 0 0 0 0 0 0 0 1 0 0 0 0
20 331 0 0 0 0 0 0 0 0 1 0 0 0
21 707 0 0 0 0 0 0 0 0 0 1 0 0
22 715 0 0 0 0 0 0 0 0 0 0 1 0
23 657 0 0 0 0 0 0 0 0 0 0 0 1
24 653 0 0 0 0 0 0 0 0 0 0 0 0
25 642 0 1 0 0 0 0 0 0 0 0 0 0
26 643 0 0 1 0 0 0 0 0 0 0 0 0
27 718 0 0 0 1 0 0 0 0 0 0 0 0
28 654 0 0 0 0 1 0 0 0 0 0 0 0
29 632 0 0 0 0 0 1 0 0 0 0 0 0
30 731 0 0 0 0 0 0 1 0 0 0 0 0
31 392 1 0 0 0 0 0 0 1 0 0 0 0
32 344 1 0 0 0 0 0 0 0 1 0 0 0
33 792 1 0 0 0 0 0 0 0 0 1 0 0
34 852 1 0 0 0 0 0 0 0 0 0 1 0
35 649 1 0 0 0 0 0 0 0 0 0 0 1
36 629 1 0 0 0 0 0 0 0 0 0 0 0
37 685 1 1 0 0 0 0 0 0 0 0 0 0
38 617 1 0 1 0 0 0 0 0 0 0 0 0
39 715 1 0 0 1 0 0 0 0 0 0 0 0
40 715 1 0 0 0 1 0 0 0 0 0 0 0
41 629 1 0 0 0 0 1 0 0 0 0 0 0
42 916 1 0 0 0 0 0 1 0 0 0 0 0
43 531 1 0 0 0 0 0 0 1 0 0 0 0
44 357 1 0 0 0 0 0 0 0 1 0 0 0
45 917 1 0 0 0 0 0 0 0 0 1 0 0
46 828 1 0 0 0 0 0 0 0 0 0 1 0
47 708 1 0 0 0 0 0 0 0 0 0 0 1
48 858 1 0 0 0 0 0 0 0 0 0 0 0
49 775 1 1 0 0 0 0 0 0 0 0 0 0
50 785 1 0 1 0 0 0 0 0 0 0 0 0
51 1006 1 0 0 1 0 0 0 0 0 0 0 0
52 789 1 0 0 0 1 0 0 0 0 0 0 0
53 734 1 0 0 0 0 1 0 0 0 0 0 0
54 906 1 0 0 0 0 0 1 0 0 0 0 0
55 532 1 0 0 0 0 0 0 1 0 0 0 0
56 387 1 0 0 0 0 0 0 0 1 0 0 0
57 991 1 0 0 0 0 0 0 0 0 1 0 0
58 841 1 0 0 0 0 0 0 0 0 0 1 0
59 892 1 0 0 0 0 0 0 0 0 0 0 1
60 782 1 0 0 0 0 0 0 0 0 0 0 0
61 813 1 1 0 0 0 0 0 0 0 0 0 0
62 793 1 0 1 0 0 0 0 0 0 0 0 0
63 978 1 0 0 1 0 0 0 0 0 0 0 0
64 775 1 0 0 0 1 0 0 0 0 0 0 0
65 797 1 0 0 0 0 1 0 0 0 0 0 0
66 946 1 0 0 0 0 0 1 0 0 0 0 0
67 594 1 0 0 0 0 0 0 1 0 0 0 0
68 438 1 0 0 0 0 0 0 0 1 0 0 0
69 1022 1 0 0 0 0 0 0 0 0 1 0 0
70 868 1 0 0 0 0 0 0 0 0 0 1 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
655.028 107.620 16.162 -14.338 101.495 -15.171
M5 M6 M7 M8 M9 M10
-1.838 138.495 -251.275 -356.108 148.725 71.725
M11
-18.200
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-149.143 -36.293 6.127 35.690 147.552
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 655.028 33.840 19.357 < 2e-16 ***
X 107.620 17.577 6.123 9.04e-08 ***
M1 16.162 43.573 0.371 0.71207
M2 -14.338 43.573 -0.329 0.74332
M3 101.495 43.573 2.329 0.02341 *
M4 -15.171 43.573 -0.348 0.72898
M5 -1.838 43.573 -0.042 0.96650
M6 138.495 43.573 3.178 0.00239 **
M7 -251.275 43.553 -5.769 3.42e-07 ***
M8 -356.108 43.553 -8.176 3.50e-11 ***
M9 148.725 43.553 3.415 0.00118 **
M10 71.725 43.553 1.647 0.10509
M11 -18.200 45.473 -0.400 0.69048
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 71.9 on 57 degrees of freedom
Multiple R-squared: 0.8444, Adjusted R-squared: 0.8117
F-statistic: 25.78 on 12 and 57 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.206978812 0.413957624 0.79302119
[2,] 0.102335650 0.204671300 0.89766435
[3,] 0.194388059 0.388776119 0.80561194
[4,] 0.115790269 0.231580538 0.88420973
[5,] 0.071106063 0.142212126 0.92889394
[6,] 0.109325743 0.218651486 0.89067426
[7,] 0.065461744 0.130923488 0.93453826
[8,] 0.044913485 0.089826970 0.95508651
[9,] 0.024930801 0.049861601 0.97506920
[10,] 0.022997123 0.045994246 0.97700288
[11,] 0.012842302 0.025684604 0.98715770
[12,] 0.006591282 0.013182564 0.99340872
[13,] 0.003832852 0.007665704 0.99616715
[14,] 0.005291551 0.010583101 0.99470845
[15,] 0.003807173 0.007614347 0.99619283
[16,] 0.003185284 0.006370569 0.99681472
[17,] 0.001619907 0.003239815 0.99838009
[18,] 0.002059146 0.004118292 0.99794085
[19,] 0.004964007 0.009928013 0.99503599
[20,] 0.004757938 0.009515875 0.99524206
[21,] 0.010228913 0.020457827 0.98977109
[22,] 0.008858453 0.017716907 0.99114155
[23,] 0.017944291 0.035888581 0.98205571
[24,] 0.233218286 0.466436572 0.76678171
[25,] 0.258504389 0.517008777 0.74149561
[26,] 0.484612213 0.969224427 0.51538779
[27,] 0.550171670 0.899656661 0.44982833
[28,] 0.579179542 0.841640916 0.42082046
[29,] 0.534971526 0.930056948 0.46502847
[30,] 0.670497444 0.659005113 0.32950256
[31,] 0.603386568 0.793226863 0.39661343
[32,] 0.951941556 0.096116888 0.04805844
[33,] 0.983056479 0.033887042 0.01694352
[34,] 0.974053274 0.051893453 0.02594673
[35,] 0.953999154 0.092001692 0.04600085
[36,] 0.958695120 0.082609761 0.04130488
[37,] 0.917386141 0.165227718 0.08261386
[38,] 0.900981479 0.198037041 0.09901852
[39,] 0.827970918 0.344058163 0.17202908
> postscript(file="/var/www/html/freestat/rcomp/tmp/1qk761292927123.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/2qk761292927123.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/3qk761292927123.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/41toq1292927123.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/51toq1292927123.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 = 70
Frequency = 1
1 2 3 4 5 6
23.809761 -2.690239 5.476428 -4.856906 67.809761 60.476428
7 8 9 10 11 12
14.246348 68.079681 20.246348 -39.753652 -35.828287 20.971713
13 14 15 16 17 18
68.809761 50.309761 -73.523572 -45.856906 75.809761 -62.523572
19 20 21 22 23 24
-17.753652 32.079681 -96.753652 -11.753652 20.171713 -2.028287
25 26 27 28 29 30
-29.190239 2.309761 -38.523572 14.143094 -21.190239 -62.523572
31 32 33 34 35 36
-119.373174 -62.539841 -119.373174 17.626826 -95.447809 -133.647809
37 38 39 40 41 42
-93.809761 -131.309761 -149.143094 -32.476428 -131.809761 14.856906
43 44 45 46 47 48
19.626826 -49.539841 5.626826 -6.373174 -36.447809 95.352191
49 50 51 52 53 54
-3.809761 36.690239 141.856906 41.523572 -26.809761 4.856906
55 56 57 58 59 60
20.626826 -19.539841 79.626826 6.626826 147.552191 19.352191
61 62 63 64 65 66
34.190239 44.690239 113.856906 27.523572 36.190239 44.856906
67 68 69 70
82.626826 31.460159 110.626826 33.626826
> postscript(file="/var/www/html/freestat/rcomp/tmp/61toq1292927123.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 = 70
Frequency = 1
lag(myerror, k = 1) myerror
0 23.809761 NA
1 -2.690239 23.809761
2 5.476428 -2.690239
3 -4.856906 5.476428
4 67.809761 -4.856906
5 60.476428 67.809761
6 14.246348 60.476428
7 68.079681 14.246348
8 20.246348 68.079681
9 -39.753652 20.246348
10 -35.828287 -39.753652
11 20.971713 -35.828287
12 68.809761 20.971713
13 50.309761 68.809761
14 -73.523572 50.309761
15 -45.856906 -73.523572
16 75.809761 -45.856906
17 -62.523572 75.809761
18 -17.753652 -62.523572
19 32.079681 -17.753652
20 -96.753652 32.079681
21 -11.753652 -96.753652
22 20.171713 -11.753652
23 -2.028287 20.171713
24 -29.190239 -2.028287
25 2.309761 -29.190239
26 -38.523572 2.309761
27 14.143094 -38.523572
28 -21.190239 14.143094
29 -62.523572 -21.190239
30 -119.373174 -62.523572
31 -62.539841 -119.373174
32 -119.373174 -62.539841
33 17.626826 -119.373174
34 -95.447809 17.626826
35 -133.647809 -95.447809
36 -93.809761 -133.647809
37 -131.309761 -93.809761
38 -149.143094 -131.309761
39 -32.476428 -149.143094
40 -131.809761 -32.476428
41 14.856906 -131.809761
42 19.626826 14.856906
43 -49.539841 19.626826
44 5.626826 -49.539841
45 -6.373174 5.626826
46 -36.447809 -6.373174
47 95.352191 -36.447809
48 -3.809761 95.352191
49 36.690239 -3.809761
50 141.856906 36.690239
51 41.523572 141.856906
52 -26.809761 41.523572
53 4.856906 -26.809761
54 20.626826 4.856906
55 -19.539841 20.626826
56 79.626826 -19.539841
57 6.626826 79.626826
58 147.552191 6.626826
59 19.352191 147.552191
60 34.190239 19.352191
61 44.690239 34.190239
62 113.856906 44.690239
63 27.523572 113.856906
64 36.190239 27.523572
65 44.856906 36.190239
66 82.626826 44.856906
67 31.460159 82.626826
68 110.626826 31.460159
69 33.626826 110.626826
70 NA 33.626826
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2.690239 23.809761
[2,] 5.476428 -2.690239
[3,] -4.856906 5.476428
[4,] 67.809761 -4.856906
[5,] 60.476428 67.809761
[6,] 14.246348 60.476428
[7,] 68.079681 14.246348
[8,] 20.246348 68.079681
[9,] -39.753652 20.246348
[10,] -35.828287 -39.753652
[11,] 20.971713 -35.828287
[12,] 68.809761 20.971713
[13,] 50.309761 68.809761
[14,] -73.523572 50.309761
[15,] -45.856906 -73.523572
[16,] 75.809761 -45.856906
[17,] -62.523572 75.809761
[18,] -17.753652 -62.523572
[19,] 32.079681 -17.753652
[20,] -96.753652 32.079681
[21,] -11.753652 -96.753652
[22,] 20.171713 -11.753652
[23,] -2.028287 20.171713
[24,] -29.190239 -2.028287
[25,] 2.309761 -29.190239
[26,] -38.523572 2.309761
[27,] 14.143094 -38.523572
[28,] -21.190239 14.143094
[29,] -62.523572 -21.190239
[30,] -119.373174 -62.523572
[31,] -62.539841 -119.373174
[32,] -119.373174 -62.539841
[33,] 17.626826 -119.373174
[34,] -95.447809 17.626826
[35,] -133.647809 -95.447809
[36,] -93.809761 -133.647809
[37,] -131.309761 -93.809761
[38,] -149.143094 -131.309761
[39,] -32.476428 -149.143094
[40,] -131.809761 -32.476428
[41,] 14.856906 -131.809761
[42,] 19.626826 14.856906
[43,] -49.539841 19.626826
[44,] 5.626826 -49.539841
[45,] -6.373174 5.626826
[46,] -36.447809 -6.373174
[47,] 95.352191 -36.447809
[48,] -3.809761 95.352191
[49,] 36.690239 -3.809761
[50,] 141.856906 36.690239
[51,] 41.523572 141.856906
[52,] -26.809761 41.523572
[53,] 4.856906 -26.809761
[54,] 20.626826 4.856906
[55,] -19.539841 20.626826
[56,] 79.626826 -19.539841
[57,] 6.626826 79.626826
[58,] 147.552191 6.626826
[59,] 19.352191 147.552191
[60,] 34.190239 19.352191
[61,] 44.690239 34.190239
[62,] 113.856906 44.690239
[63,] 27.523572 113.856906
[64,] 36.190239 27.523572
[65,] 44.856906 36.190239
[66,] 82.626826 44.856906
[67,] 31.460159 82.626826
[68,] 110.626826 31.460159
[69,] 33.626826 110.626826
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2.690239 23.809761
2 5.476428 -2.690239
3 -4.856906 5.476428
4 67.809761 -4.856906
5 60.476428 67.809761
6 14.246348 60.476428
7 68.079681 14.246348
8 20.246348 68.079681
9 -39.753652 20.246348
10 -35.828287 -39.753652
11 20.971713 -35.828287
12 68.809761 20.971713
13 50.309761 68.809761
14 -73.523572 50.309761
15 -45.856906 -73.523572
16 75.809761 -45.856906
17 -62.523572 75.809761
18 -17.753652 -62.523572
19 32.079681 -17.753652
20 -96.753652 32.079681
21 -11.753652 -96.753652
22 20.171713 -11.753652
23 -2.028287 20.171713
24 -29.190239 -2.028287
25 2.309761 -29.190239
26 -38.523572 2.309761
27 14.143094 -38.523572
28 -21.190239 14.143094
29 -62.523572 -21.190239
30 -119.373174 -62.523572
31 -62.539841 -119.373174
32 -119.373174 -62.539841
33 17.626826 -119.373174
34 -95.447809 17.626826
35 -133.647809 -95.447809
36 -93.809761 -133.647809
37 -131.309761 -93.809761
38 -149.143094 -131.309761
39 -32.476428 -149.143094
40 -131.809761 -32.476428
41 14.856906 -131.809761
42 19.626826 14.856906
43 -49.539841 19.626826
44 5.626826 -49.539841
45 -6.373174 5.626826
46 -36.447809 -6.373174
47 95.352191 -36.447809
48 -3.809761 95.352191
49 36.690239 -3.809761
50 141.856906 36.690239
51 41.523572 141.856906
52 -26.809761 41.523572
53 4.856906 -26.809761
54 20.626826 4.856906
55 -19.539841 20.626826
56 79.626826 -19.539841
57 6.626826 79.626826
58 147.552191 6.626826
59 19.352191 147.552191
60 34.190239 19.352191
61 44.690239 34.190239
62 113.856906 44.690239
63 27.523572 113.856906
64 36.190239 27.523572
65 44.856906 36.190239
66 82.626826 44.856906
67 31.460159 82.626826
68 110.626826 31.460159
69 33.626826 110.626826
> 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/7u3ou1292927123.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/84une1292927123.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/94une1292927123.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/104une1292927123.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/110l251292927123.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/12bvk81292927123.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/13iwz21292927123.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/14snyn1292927123.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/15wofb1292927123.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/16sfuj1292927123.tab")
+ }
>
> try(system("convert tmp/1qk761292927123.ps tmp/1qk761292927123.png",intern=TRUE))
character(0)
> try(system("convert tmp/2qk761292927123.ps tmp/2qk761292927123.png",intern=TRUE))
character(0)
> try(system("convert tmp/3qk761292927123.ps tmp/3qk761292927123.png",intern=TRUE))
character(0)
> try(system("convert tmp/41toq1292927123.ps tmp/41toq1292927123.png",intern=TRUE))
character(0)
> try(system("convert tmp/51toq1292927123.ps tmp/51toq1292927123.png",intern=TRUE))
character(0)
> try(system("convert tmp/61toq1292927123.ps tmp/61toq1292927123.png",intern=TRUE))
character(0)
> try(system("convert tmp/7u3ou1292927123.ps tmp/7u3ou1292927123.png",intern=TRUE))
character(0)
> try(system("convert tmp/84une1292927123.ps tmp/84une1292927123.png",intern=TRUE))
character(0)
> try(system("convert tmp/94une1292927123.ps tmp/94une1292927123.png",intern=TRUE))
character(0)
> try(system("convert tmp/104une1292927123.ps tmp/104une1292927123.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.032 2.582 4.342