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(267413,21.4,267366,26.4,264777,26.4,258863,29.4,254844,34.4,254868,24.4,277267,26.4,285351,25.4,286602,31.4,283042,27.4,276687,27.4,277915,29.4,277128,32.4,277103,26.4,275037,22.4,270150,19.4,267140,21.4,264993,23.4,287259,23.4,291186,25.4,292300,28.4,288186,27.4,281477,21.4,282656,17.4,280190,24.4,280408,26.4,276836,22.4,275216,14.4,274352,18.4,271311,25.4,289802,29.4,290726,26.4,292300,26.4,278506,20.4,269826,26.4,265861,29.4,269034,33.4,264176,32.4,255198,35.4,253353,34.4,246057,36.4,235372,32.4,258556,34.4,260993,31.4,254663,27.4,250643,27.4,243422,30.4,247105,32.4,248541,32.4,245039,27.4,237080,31.4,237085,29.4,225554,27.4,226839,25.4,247934,26.4,248333,23.4,246969,18.4,245098,22.4,246263,17.4,255765,17.4,264319,11.4,268347,9.4,273046,6.4,273963,0,267430,7.8,271993,7.9,292710,12,295881,16.9,293299,12.3),dim=c(2,69),dimnames=list(c('Y','X'),1:69))
> y <- array(NA,dim=c(2,69),dimnames=list(c('Y','X'),1:69))
> 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 267413 21.4 1 0 0 0 0 0 0 0 0 0 0
2 267366 26.4 0 1 0 0 0 0 0 0 0 0 0
3 264777 26.4 0 0 1 0 0 0 0 0 0 0 0
4 258863 29.4 0 0 0 1 0 0 0 0 0 0 0
5 254844 34.4 0 0 0 0 1 0 0 0 0 0 0
6 254868 24.4 0 0 0 0 0 1 0 0 0 0 0
7 277267 26.4 0 0 0 0 0 0 1 0 0 0 0
8 285351 25.4 0 0 0 0 0 0 0 1 0 0 0
9 286602 31.4 0 0 0 0 0 0 0 0 1 0 0
10 283042 27.4 0 0 0 0 0 0 0 0 0 1 0
11 276687 27.4 0 0 0 0 0 0 0 0 0 0 1
12 277915 29.4 0 0 0 0 0 0 0 0 0 0 0
13 277128 32.4 1 0 0 0 0 0 0 0 0 0 0
14 277103 26.4 0 1 0 0 0 0 0 0 0 0 0
15 275037 22.4 0 0 1 0 0 0 0 0 0 0 0
16 270150 19.4 0 0 0 1 0 0 0 0 0 0 0
17 267140 21.4 0 0 0 0 1 0 0 0 0 0 0
18 264993 23.4 0 0 0 0 0 1 0 0 0 0 0
19 287259 23.4 0 0 0 0 0 0 1 0 0 0 0
20 291186 25.4 0 0 0 0 0 0 0 1 0 0 0
21 292300 28.4 0 0 0 0 0 0 0 0 1 0 0
22 288186 27.4 0 0 0 0 0 0 0 0 0 1 0
23 281477 21.4 0 0 0 0 0 0 0 0 0 0 1
24 282656 17.4 0 0 0 0 0 0 0 0 0 0 0
25 280190 24.4 1 0 0 0 0 0 0 0 0 0 0
26 280408 26.4 0 1 0 0 0 0 0 0 0 0 0
27 276836 22.4 0 0 1 0 0 0 0 0 0 0 0
28 275216 14.4 0 0 0 1 0 0 0 0 0 0 0
29 274352 18.4 0 0 0 0 1 0 0 0 0 0 0
30 271311 25.4 0 0 0 0 0 1 0 0 0 0 0
31 289802 29.4 0 0 0 0 0 0 1 0 0 0 0
32 290726 26.4 0 0 0 0 0 0 0 1 0 0 0
33 292300 26.4 0 0 0 0 0 0 0 0 1 0 0
34 278506 20.4 0 0 0 0 0 0 0 0 0 1 0
35 269826 26.4 0 0 0 0 0 0 0 0 0 0 1
36 265861 29.4 0 0 0 0 0 0 0 0 0 0 0
37 269034 33.4 1 0 0 0 0 0 0 0 0 0 0
38 264176 32.4 0 1 0 0 0 0 0 0 0 0 0
39 255198 35.4 0 0 1 0 0 0 0 0 0 0 0
40 253353 34.4 0 0 0 1 0 0 0 0 0 0 0
41 246057 36.4 0 0 0 0 1 0 0 0 0 0 0
42 235372 32.4 0 0 0 0 0 1 0 0 0 0 0
43 258556 34.4 0 0 0 0 0 0 1 0 0 0 0
44 260993 31.4 0 0 0 0 0 0 0 1 0 0 0
45 254663 27.4 0 0 0 0 0 0 0 0 1 0 0
46 250643 27.4 0 0 0 0 0 0 0 0 0 1 0
47 243422 30.4 0 0 0 0 0 0 0 0 0 0 1
48 247105 32.4 0 0 0 0 0 0 0 0 0 0 0
49 248541 32.4 1 0 0 0 0 0 0 0 0 0 0
50 245039 27.4 0 1 0 0 0 0 0 0 0 0 0
51 237080 31.4 0 0 1 0 0 0 0 0 0 0 0
52 237085 29.4 0 0 0 1 0 0 0 0 0 0 0
53 225554 27.4 0 0 0 0 1 0 0 0 0 0 0
54 226839 25.4 0 0 0 0 0 1 0 0 0 0 0
55 247934 26.4 0 0 0 0 0 0 1 0 0 0 0
56 248333 23.4 0 0 0 0 0 0 0 1 0 0 0
57 246969 18.4 0 0 0 0 0 0 0 0 1 0 0
58 245098 22.4 0 0 0 0 0 0 0 0 0 1 0
59 246263 17.4 0 0 0 0 0 0 0 0 0 0 1
60 255765 17.4 0 0 0 0 0 0 0 0 0 0 0
61 264319 11.4 1 0 0 0 0 0 0 0 0 0 0
62 268347 9.4 0 1 0 0 0 0 0 0 0 0 0
63 273046 6.4 0 0 1 0 0 0 0 0 0 0 0
64 273963 0.0 0 0 0 1 0 0 0 0 0 0 0
65 267430 7.8 0 0 0 0 1 0 0 0 0 0 0
66 271993 7.9 0 0 0 0 0 1 0 0 0 0 0
67 292710 12.0 0 0 0 0 0 0 1 0 0 0 0
68 295881 16.9 0 0 0 0 0 0 0 1 0 0 0
69 293299 12.3 0 0 0 0 0 0 0 0 1 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
284296.2 -731.6 2422.5 871.4 -3027.2 -7372.8
M5 M6 M7 M8 M9 M10
-10622.7 -13130.8 9825.1 12604.2 10987.1 3088.3
M11
-2764.3
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-34853 -12936 3448 11322 20847
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 284296.2 9679.2 29.372 < 2e-16 ***
X -731.6 257.1 -2.846 0.00618 **
M1 2422.5 9739.3 0.249 0.80447
M2 871.4 9738.4 0.089 0.92902
M3 -3027.2 9742.0 -0.311 0.75716
M4 -7372.8 9792.7 -0.753 0.45467
M5 -10622.7 9740.4 -1.091 0.28013
M6 -13130.8 9751.9 -1.346 0.18357
M7 9825.1 9737.7 1.009 0.31733
M8 12604.2 9738.2 1.294 0.20087
M9 10987.1 9742.2 1.128 0.26422
M10 3088.3 10170.8 0.304 0.76253
M11 -2764.3 10171.8 -0.272 0.78680
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 16080 on 56 degrees of freedom
Multiple R-squared: 0.2967, Adjusted R-squared: 0.1459
F-statistic: 1.968 on 12 and 56 DF, p-value: 0.04506
> 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,] 1.182242e-01 2.364484e-01 0.8817758
[2,] 5.492340e-02 1.098468e-01 0.9450766
[3,] 2.947866e-02 5.895732e-02 0.9705213
[4,] 1.449413e-02 2.898826e-02 0.9855059
[5,] 6.287626e-03 1.257525e-02 0.9937124
[6,] 2.709692e-03 5.419385e-03 0.9972903
[7,] 1.395030e-03 2.790060e-03 0.9986050
[8,] 5.533276e-04 1.106655e-03 0.9994467
[9,] 2.020949e-04 4.041898e-04 0.9997979
[10,] 1.029460e-04 2.058920e-04 0.9998971
[11,] 7.530102e-05 1.506020e-04 0.9999247
[12,] 3.797034e-05 7.594067e-05 0.9999620
[13,] 1.489539e-05 2.979078e-05 0.9999851
[14,] 8.840121e-06 1.768024e-05 0.9999912
[15,] 2.008186e-05 4.016373e-05 0.9999799
[16,] 3.044264e-05 6.088528e-05 0.9999696
[17,] 1.900520e-05 3.801040e-05 0.9999810
[18,] 2.080318e-05 4.160636e-05 0.9999792
[19,] 4.082593e-05 8.165186e-05 0.9999592
[20,] 5.574164e-05 1.114833e-04 0.9999443
[21,] 6.833639e-05 1.366728e-04 0.9999317
[22,] 6.040245e-05 1.208049e-04 0.9999396
[23,] 6.995755e-05 1.399151e-04 0.9999300
[24,] 7.728825e-05 1.545765e-04 0.9999227
[25,] 8.354343e-05 1.670869e-04 0.9999165
[26,] 2.041744e-04 4.083489e-04 0.9997958
[27,] 8.548518e-04 1.709704e-03 0.9991451
[28,] 1.608089e-03 3.216179e-03 0.9983919
[29,] 3.374045e-03 6.748091e-03 0.9966260
[30,] 2.667099e-02 5.334197e-02 0.9733290
[31,] 4.978829e-02 9.957658e-02 0.9502117
[32,] 7.034969e-02 1.406994e-01 0.9296503
[33,] 7.009478e-02 1.401896e-01 0.9299052
[34,] 8.425902e-02 1.685180e-01 0.9157410
[35,] 8.767949e-02 1.753590e-01 0.9123205
[36,] 8.124855e-02 1.624971e-01 0.9187515
[37,] 2.008540e-01 4.017079e-01 0.7991460
[38,] 2.373528e-01 4.747057e-01 0.7626472
> postscript(file="/var/www/html/rcomp/tmp/1fwjr1258714437.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/2jwow1258714437.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/3tzhw1258714437.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/4z6rd1258714437.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/5zady1258714437.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 = 69
Frequency = 1
1 2 3 4 5 6
-3649.9424 1512.1330 2821.6862 3448.0069 6336.7892 1553.1414
7 8 9 10 11 12
2459.3518 7032.7549 14290.2781 15702.7915 15200.4234 15127.2351
13 14 15 16 17 18
14112.4353 11249.1330 10155.3670 7419.2090 9122.2519 10946.5616
19 20 21 22 23 24
10256.6124 12867.7549 17793.5387 20846.7915 15600.9447 11089.2776
25 26 27 28 29 30
11321.7970 14554.1330 11954.3670 8827.3101 14139.5126 18727.7212
31 32 33 34 35 36
17189.0911 13139.3347 16330.3792 6045.7330 7607.8436 3073.2351
37 38 39 40 41 42
6750.0151 2711.6117 -173.0957 1595.9059 -987.0512 -12090.2203
43 44 45 46 47 48
-10399.0099 -12935.7664 -20575.0410 -16696.2085 -15869.8372 -13488.0255
49 50 51 52 53 54
-14474.5647 -20083.2872 -21217.4149 -18329.9931 -28074.2693 -25744.2788
55 56 57 58 59 60
-26873.6482 -31448.4047 -34853.2591 -25899.1075 -22539.3745 -15801.7224
61 62 63 64 65 66
-14059.7403 -9943.7234 -3540.9096 -2960.4389 -537.2332 6607.0749
67 68 69
7367.6028 11344.3267 7014.1041
> postscript(file="/var/www/html/rcomp/tmp/6fjrl1258714437.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 = 69
Frequency = 1
lag(myerror, k = 1) myerror
0 -3649.9424 NA
1 1512.1330 -3649.9424
2 2821.6862 1512.1330
3 3448.0069 2821.6862
4 6336.7892 3448.0069
5 1553.1414 6336.7892
6 2459.3518 1553.1414
7 7032.7549 2459.3518
8 14290.2781 7032.7549
9 15702.7915 14290.2781
10 15200.4234 15702.7915
11 15127.2351 15200.4234
12 14112.4353 15127.2351
13 11249.1330 14112.4353
14 10155.3670 11249.1330
15 7419.2090 10155.3670
16 9122.2519 7419.2090
17 10946.5616 9122.2519
18 10256.6124 10946.5616
19 12867.7549 10256.6124
20 17793.5387 12867.7549
21 20846.7915 17793.5387
22 15600.9447 20846.7915
23 11089.2776 15600.9447
24 11321.7970 11089.2776
25 14554.1330 11321.7970
26 11954.3670 14554.1330
27 8827.3101 11954.3670
28 14139.5126 8827.3101
29 18727.7212 14139.5126
30 17189.0911 18727.7212
31 13139.3347 17189.0911
32 16330.3792 13139.3347
33 6045.7330 16330.3792
34 7607.8436 6045.7330
35 3073.2351 7607.8436
36 6750.0151 3073.2351
37 2711.6117 6750.0151
38 -173.0957 2711.6117
39 1595.9059 -173.0957
40 -987.0512 1595.9059
41 -12090.2203 -987.0512
42 -10399.0099 -12090.2203
43 -12935.7664 -10399.0099
44 -20575.0410 -12935.7664
45 -16696.2085 -20575.0410
46 -15869.8372 -16696.2085
47 -13488.0255 -15869.8372
48 -14474.5647 -13488.0255
49 -20083.2872 -14474.5647
50 -21217.4149 -20083.2872
51 -18329.9931 -21217.4149
52 -28074.2693 -18329.9931
53 -25744.2788 -28074.2693
54 -26873.6482 -25744.2788
55 -31448.4047 -26873.6482
56 -34853.2591 -31448.4047
57 -25899.1075 -34853.2591
58 -22539.3745 -25899.1075
59 -15801.7224 -22539.3745
60 -14059.7403 -15801.7224
61 -9943.7234 -14059.7403
62 -3540.9096 -9943.7234
63 -2960.4389 -3540.9096
64 -537.2332 -2960.4389
65 6607.0749 -537.2332
66 7367.6028 6607.0749
67 11344.3267 7367.6028
68 7014.1041 11344.3267
69 NA 7014.1041
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1512.1330 -3649.9424
[2,] 2821.6862 1512.1330
[3,] 3448.0069 2821.6862
[4,] 6336.7892 3448.0069
[5,] 1553.1414 6336.7892
[6,] 2459.3518 1553.1414
[7,] 7032.7549 2459.3518
[8,] 14290.2781 7032.7549
[9,] 15702.7915 14290.2781
[10,] 15200.4234 15702.7915
[11,] 15127.2351 15200.4234
[12,] 14112.4353 15127.2351
[13,] 11249.1330 14112.4353
[14,] 10155.3670 11249.1330
[15,] 7419.2090 10155.3670
[16,] 9122.2519 7419.2090
[17,] 10946.5616 9122.2519
[18,] 10256.6124 10946.5616
[19,] 12867.7549 10256.6124
[20,] 17793.5387 12867.7549
[21,] 20846.7915 17793.5387
[22,] 15600.9447 20846.7915
[23,] 11089.2776 15600.9447
[24,] 11321.7970 11089.2776
[25,] 14554.1330 11321.7970
[26,] 11954.3670 14554.1330
[27,] 8827.3101 11954.3670
[28,] 14139.5126 8827.3101
[29,] 18727.7212 14139.5126
[30,] 17189.0911 18727.7212
[31,] 13139.3347 17189.0911
[32,] 16330.3792 13139.3347
[33,] 6045.7330 16330.3792
[34,] 7607.8436 6045.7330
[35,] 3073.2351 7607.8436
[36,] 6750.0151 3073.2351
[37,] 2711.6117 6750.0151
[38,] -173.0957 2711.6117
[39,] 1595.9059 -173.0957
[40,] -987.0512 1595.9059
[41,] -12090.2203 -987.0512
[42,] -10399.0099 -12090.2203
[43,] -12935.7664 -10399.0099
[44,] -20575.0410 -12935.7664
[45,] -16696.2085 -20575.0410
[46,] -15869.8372 -16696.2085
[47,] -13488.0255 -15869.8372
[48,] -14474.5647 -13488.0255
[49,] -20083.2872 -14474.5647
[50,] -21217.4149 -20083.2872
[51,] -18329.9931 -21217.4149
[52,] -28074.2693 -18329.9931
[53,] -25744.2788 -28074.2693
[54,] -26873.6482 -25744.2788
[55,] -31448.4047 -26873.6482
[56,] -34853.2591 -31448.4047
[57,] -25899.1075 -34853.2591
[58,] -22539.3745 -25899.1075
[59,] -15801.7224 -22539.3745
[60,] -14059.7403 -15801.7224
[61,] -9943.7234 -14059.7403
[62,] -3540.9096 -9943.7234
[63,] -2960.4389 -3540.9096
[64,] -537.2332 -2960.4389
[65,] 6607.0749 -537.2332
[66,] 7367.6028 6607.0749
[67,] 11344.3267 7367.6028
[68,] 7014.1041 11344.3267
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1512.1330 -3649.9424
2 2821.6862 1512.1330
3 3448.0069 2821.6862
4 6336.7892 3448.0069
5 1553.1414 6336.7892
6 2459.3518 1553.1414
7 7032.7549 2459.3518
8 14290.2781 7032.7549
9 15702.7915 14290.2781
10 15200.4234 15702.7915
11 15127.2351 15200.4234
12 14112.4353 15127.2351
13 11249.1330 14112.4353
14 10155.3670 11249.1330
15 7419.2090 10155.3670
16 9122.2519 7419.2090
17 10946.5616 9122.2519
18 10256.6124 10946.5616
19 12867.7549 10256.6124
20 17793.5387 12867.7549
21 20846.7915 17793.5387
22 15600.9447 20846.7915
23 11089.2776 15600.9447
24 11321.7970 11089.2776
25 14554.1330 11321.7970
26 11954.3670 14554.1330
27 8827.3101 11954.3670
28 14139.5126 8827.3101
29 18727.7212 14139.5126
30 17189.0911 18727.7212
31 13139.3347 17189.0911
32 16330.3792 13139.3347
33 6045.7330 16330.3792
34 7607.8436 6045.7330
35 3073.2351 7607.8436
36 6750.0151 3073.2351
37 2711.6117 6750.0151
38 -173.0957 2711.6117
39 1595.9059 -173.0957
40 -987.0512 1595.9059
41 -12090.2203 -987.0512
42 -10399.0099 -12090.2203
43 -12935.7664 -10399.0099
44 -20575.0410 -12935.7664
45 -16696.2085 -20575.0410
46 -15869.8372 -16696.2085
47 -13488.0255 -15869.8372
48 -14474.5647 -13488.0255
49 -20083.2872 -14474.5647
50 -21217.4149 -20083.2872
51 -18329.9931 -21217.4149
52 -28074.2693 -18329.9931
53 -25744.2788 -28074.2693
54 -26873.6482 -25744.2788
55 -31448.4047 -26873.6482
56 -34853.2591 -31448.4047
57 -25899.1075 -34853.2591
58 -22539.3745 -25899.1075
59 -15801.7224 -22539.3745
60 -14059.7403 -15801.7224
61 -9943.7234 -14059.7403
62 -3540.9096 -9943.7234
63 -2960.4389 -3540.9096
64 -537.2332 -2960.4389
65 6607.0749 -537.2332
66 7367.6028 6607.0749
67 11344.3267 7367.6028
68 7014.1041 11344.3267
> 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/7q8u41258714437.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/86izn1258714437.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/9tyjn1258714437.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/10ix5o1258714437.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/113s4a1258714437.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/123bah1258714437.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/13368i1258714437.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/142sbh1258714437.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/15dp2v1258714437.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/16cbuu1258714437.tab")
+ }
>
> system("convert tmp/1fwjr1258714437.ps tmp/1fwjr1258714437.png")
> system("convert tmp/2jwow1258714437.ps tmp/2jwow1258714437.png")
> system("convert tmp/3tzhw1258714437.ps tmp/3tzhw1258714437.png")
> system("convert tmp/4z6rd1258714437.ps tmp/4z6rd1258714437.png")
> system("convert tmp/5zady1258714437.ps tmp/5zady1258714437.png")
> system("convert tmp/6fjrl1258714437.ps tmp/6fjrl1258714437.png")
> system("convert tmp/7q8u41258714437.ps tmp/7q8u41258714437.png")
> system("convert tmp/86izn1258714437.ps tmp/86izn1258714437.png")
> system("convert tmp/9tyjn1258714437.ps tmp/9tyjn1258714437.png")
> system("convert tmp/10ix5o1258714437.ps tmp/10ix5o1258714437.png")
>
>
> proc.time()
user system elapsed
2.469 1.561 2.983