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(3884.3
+ ,12476.8
+ ,3983.4
+ ,3956.2
+ ,3892.2
+ ,12384.6
+ ,4152.9
+ ,3142.7
+ ,3613
+ ,12266.7
+ ,4286.1
+ ,3884.3
+ ,3730.5
+ ,12919.9
+ ,4348.1
+ ,3892.2
+ ,3481.3
+ ,11497.3
+ ,3949.3
+ ,3613
+ ,3649.5
+ ,12142
+ ,4166.7
+ ,3730.5
+ ,4215.2
+ ,13919.4
+ ,4217.9
+ ,3481.3
+ ,4066.6
+ ,12656.8
+ ,4528.2
+ ,3649.5
+ ,4196.8
+ ,12034.1
+ ,4232.2
+ ,4215.2
+ ,4536.6
+ ,13199.7
+ ,4470.9
+ ,4066.6
+ ,4441.6
+ ,10881.3
+ ,5121.2
+ ,4196.8
+ ,3548.3
+ ,11301.2
+ ,4170.8
+ ,4536.6
+ ,4735.9
+ ,13643.9
+ ,4398.6
+ ,4441.6
+ ,4130.6
+ ,12517
+ ,4491.4
+ ,3548.3
+ ,4356.2
+ ,13981.1
+ ,4251.8
+ ,4735.9
+ ,4159.6
+ ,14275.7
+ ,4901.9
+ ,4130.6
+ ,3988
+ ,13435
+ ,4745.2
+ ,4356.2
+ ,4167.8
+ ,13565.7
+ ,4666.9
+ ,4159.6
+ ,4902.2
+ ,16216.3
+ ,4210.4
+ ,3988
+ ,3909.4
+ ,12970
+ ,5273.6
+ ,4167.8
+ ,4697.6
+ ,14079.9
+ ,4095.3
+ ,4902.2
+ ,4308.9
+ ,14235
+ ,4610.1
+ ,3909.4
+ ,4420.4
+ ,12213.4
+ ,4718.1
+ ,4697.6
+ ,3544.2
+ ,12581
+ ,4185.5
+ ,4308.9
+ ,4433
+ ,14130.4
+ ,4314.7
+ ,4420.4
+ ,4479.7
+ ,14210.8
+ ,4422.6
+ ,3544.2
+ ,4533.2
+ ,14378.5
+ ,5059.2
+ ,4433
+ ,4237.5
+ ,13142.8
+ ,5043.6
+ ,4479.7
+ ,4207.4
+ ,13714.7
+ ,4436.6
+ ,4533.2
+ ,4394
+ ,13621.9
+ ,4922.6
+ ,4237.5
+ ,5148.4
+ ,15379.8
+ ,4454.8
+ ,4207.4
+ ,4202.2
+ ,13306.3
+ ,5058.7
+ ,4394
+ ,4682.5
+ ,14391.2
+ ,4768.9
+ ,5148.4
+ ,4884.3
+ ,14909.9
+ ,5171.8
+ ,4202.2
+ ,5288.9
+ ,14025.4
+ ,4989.3
+ ,4682.5
+ ,4505.2
+ ,12951.2
+ ,5202.1
+ ,4884.3
+ ,4611.5
+ ,14344.3
+ ,4838.4
+ ,5288.9
+ ,5104
+ ,16093.4
+ ,4876.5
+ ,4505.2
+ ,4586.6
+ ,15413.6
+ ,5875.5
+ ,4611.5
+ ,4529.3
+ ,14705.7
+ ,5717.9
+ ,5104
+ ,4504.1
+ ,15972.8
+ ,4778.8
+ ,4586.6
+ ,4604.9
+ ,16241.4
+ ,6195.9
+ ,4529.3
+ ,4795.4
+ ,16626.4
+ ,4625.4
+ ,4504.1
+ ,5391.1
+ ,17136.2
+ ,5549.8
+ ,4604.9
+ ,5213.9
+ ,15622.9
+ ,6397.6
+ ,4795.4
+ ,5415
+ ,18003.9
+ ,5856.7
+ ,5391.1
+ ,5990.3
+ ,16136.1
+ ,6343.8
+ ,5213.9
+ ,4241.8
+ ,14423.7
+ ,6615.5
+ ,5415
+ ,5677.6
+ ,16789.4
+ ,5904.6
+ ,5990.3
+ ,5164.2
+ ,16782.2
+ ,6861
+ ,4241.8
+ ,3962.3
+ ,14133.8
+ ,6553.5
+ ,5677.6
+ ,4011
+ ,12607
+ ,5481
+ ,5164.2
+ ,3310.3
+ ,12004.5
+ ,5435.3
+ ,3962.3
+ ,3837.3
+ ,12175.4
+ ,5278
+ ,4011
+ ,4145.3
+ ,13268
+ ,4671.8
+ ,3310.3
+ ,3796.7
+ ,12299.3
+ ,4891.5
+ ,3837.3
+ ,3849.6
+ ,11800.6
+ ,4241.6
+ ,4145.3
+ ,4285
+ ,13873.3
+ ,4152.1
+ ,3796.7
+ ,4189.6
+ ,12269.6
+ ,4484.4
+ ,3849.6)
+ ,dim=c(4
+ ,59)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y1'
+ ,'Y2')
+ ,1:59))
> y <- array(NA,dim=c(4,59),dimnames=list(c('Y','X','Y1','Y2'),1:59))
> 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 Y1 Y2 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 3884.3 12476.8 3983.4 3956.2 1 0 0 0 0 0 0 0 0 0 0 1
2 3892.2 12384.6 4152.9 3142.7 0 1 0 0 0 0 0 0 0 0 0 2
3 3613.0 12266.7 4286.1 3884.3 0 0 1 0 0 0 0 0 0 0 0 3
4 3730.5 12919.9 4348.1 3892.2 0 0 0 1 0 0 0 0 0 0 0 4
5 3481.3 11497.3 3949.3 3613.0 0 0 0 0 1 0 0 0 0 0 0 5
6 3649.5 12142.0 4166.7 3730.5 0 0 0 0 0 1 0 0 0 0 0 6
7 4215.2 13919.4 4217.9 3481.3 0 0 0 0 0 0 1 0 0 0 0 7
8 4066.6 12656.8 4528.2 3649.5 0 0 0 0 0 0 0 1 0 0 0 8
9 4196.8 12034.1 4232.2 4215.2 0 0 0 0 0 0 0 0 1 0 0 9
10 4536.6 13199.7 4470.9 4066.6 0 0 0 0 0 0 0 0 0 1 0 10
11 4441.6 10881.3 5121.2 4196.8 0 0 0 0 0 0 0 0 0 0 1 11
12 3548.3 11301.2 4170.8 4536.6 0 0 0 0 0 0 0 0 0 0 0 12
13 4735.9 13643.9 4398.6 4441.6 1 0 0 0 0 0 0 0 0 0 0 13
14 4130.6 12517.0 4491.4 3548.3 0 1 0 0 0 0 0 0 0 0 0 14
15 4356.2 13981.1 4251.8 4735.9 0 0 1 0 0 0 0 0 0 0 0 15
16 4159.6 14275.7 4901.9 4130.6 0 0 0 1 0 0 0 0 0 0 0 16
17 3988.0 13435.0 4745.2 4356.2 0 0 0 0 1 0 0 0 0 0 0 17
18 4167.8 13565.7 4666.9 4159.6 0 0 0 0 0 1 0 0 0 0 0 18
19 4902.2 16216.3 4210.4 3988.0 0 0 0 0 0 0 1 0 0 0 0 19
20 3909.4 12970.0 5273.6 4167.8 0 0 0 0 0 0 0 1 0 0 0 20
21 4697.6 14079.9 4095.3 4902.2 0 0 0 0 0 0 0 0 1 0 0 21
22 4308.9 14235.0 4610.1 3909.4 0 0 0 0 0 0 0 0 0 1 0 22
23 4420.4 12213.4 4718.1 4697.6 0 0 0 0 0 0 0 0 0 0 1 23
24 3544.2 12581.0 4185.5 4308.9 0 0 0 0 0 0 0 0 0 0 0 24
25 4433.0 14130.4 4314.7 4420.4 1 0 0 0 0 0 0 0 0 0 0 25
26 4479.7 14210.8 4422.6 3544.2 0 1 0 0 0 0 0 0 0 0 0 26
27 4533.2 14378.5 5059.2 4433.0 0 0 1 0 0 0 0 0 0 0 0 27
28 4237.5 13142.8 5043.6 4479.7 0 0 0 1 0 0 0 0 0 0 0 28
29 4207.4 13714.7 4436.6 4533.2 0 0 0 0 1 0 0 0 0 0 0 29
30 4394.0 13621.9 4922.6 4237.5 0 0 0 0 0 1 0 0 0 0 0 30
31 5148.4 15379.8 4454.8 4207.4 0 0 0 0 0 0 1 0 0 0 0 31
32 4202.2 13306.3 5058.7 4394.0 0 0 0 0 0 0 0 1 0 0 0 32
33 4682.5 14391.2 4768.9 5148.4 0 0 0 0 0 0 0 0 1 0 0 33
34 4884.3 14909.9 5171.8 4202.2 0 0 0 0 0 0 0 0 0 1 0 34
35 5288.9 14025.4 4989.3 4682.5 0 0 0 0 0 0 0 0 0 0 1 35
36 4505.2 12951.2 5202.1 4884.3 0 0 0 0 0 0 0 0 0 0 0 36
37 4611.5 14344.3 4838.4 5288.9 1 0 0 0 0 0 0 0 0 0 0 37
38 5104.0 16093.4 4876.5 4505.2 0 1 0 0 0 0 0 0 0 0 0 38
39 4586.6 15413.6 5875.5 4611.5 0 0 1 0 0 0 0 0 0 0 0 39
40 4529.3 14705.7 5717.9 5104.0 0 0 0 1 0 0 0 0 0 0 0 40
41 4504.1 15972.8 4778.8 4586.6 0 0 0 0 1 0 0 0 0 0 0 41
42 4604.9 16241.4 6195.9 4529.3 0 0 0 0 0 1 0 0 0 0 0 42
43 4795.4 16626.4 4625.4 4504.1 0 0 0 0 0 0 1 0 0 0 0 43
44 5391.1 17136.2 5549.8 4604.9 0 0 0 0 0 0 0 1 0 0 0 44
45 5213.9 15622.9 6397.6 4795.4 0 0 0 0 0 0 0 0 1 0 0 45
46 5415.0 18003.9 5856.7 5391.1 0 0 0 0 0 0 0 0 0 1 0 46
47 5990.3 16136.1 6343.8 5213.9 0 0 0 0 0 0 0 0 0 0 1 47
48 4241.8 14423.7 6615.5 5415.0 0 0 0 0 0 0 0 0 0 0 0 48
49 5677.6 16789.4 5904.6 5990.3 1 0 0 0 0 0 0 0 0 0 0 49
50 5164.2 16782.2 6861.0 4241.8 0 1 0 0 0 0 0 0 0 0 0 50
51 3962.3 14133.8 6553.5 5677.6 0 0 1 0 0 0 0 0 0 0 0 51
52 4011.0 12607.0 5481.0 5164.2 0 0 0 1 0 0 0 0 0 0 0 52
53 3310.3 12004.5 5435.3 3962.3 0 0 0 0 1 0 0 0 0 0 0 53
54 3837.3 12175.4 5278.0 4011.0 0 0 0 0 0 1 0 0 0 0 0 54
55 4145.3 13268.0 4671.8 3310.3 0 0 0 0 0 0 1 0 0 0 0 55
56 3796.7 12299.3 4891.5 3837.3 0 0 0 0 0 0 0 1 0 0 0 56
57 3849.6 11800.6 4241.6 4145.3 0 0 0 0 0 0 0 0 1 0 0 57
58 4285.0 13873.3 4152.1 3796.7 0 0 0 0 0 0 0 0 0 1 0 58
59 4189.6 12269.6 4484.4 3849.6 0 0 0 0 0 0 0 0 0 0 1 59
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X Y1 Y2 M1 M2
-169.54897 0.25580 0.02073 0.17589 320.15062 352.41918
M3 M4 M5 M6 M7 M8
-53.92813 23.92602 -86.35697 97.55154 273.29281 214.74977
M9 M10 M11 t
411.97643 313.63195 891.25787 -3.16007
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-402.90 -125.36 -1.74 141.85 508.72
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -169.54897 331.38634 -0.512 0.61152
X 0.25580 0.03084 8.294 1.82e-10 ***
Y1 0.02073 0.06703 0.309 0.75865
Y2 0.17589 0.09889 1.779 0.08235 .
M1 320.15062 156.23455 2.049 0.04658 *
M2 352.41918 195.81067 1.800 0.07891 .
M3 -53.92813 153.44050 -0.351 0.72696
M4 23.92602 151.85553 0.158 0.87554
M5 -86.35697 160.75380 -0.537 0.59390
M6 97.55154 166.59948 0.586 0.56124
M7 273.29281 204.91515 1.334 0.18933
M8 214.74977 168.44607 1.275 0.20920
M9 411.97643 152.58014 2.700 0.00987 **
M10 313.63195 178.80378 1.754 0.08655 .
M11 891.25787 150.59395 5.918 4.81e-07 ***
t -3.16007 2.20560 -1.433 0.15916
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 218.4 on 43 degrees of freedom
Multiple R-squared: 0.8885, Adjusted R-squared: 0.8496
F-statistic: 22.85 on 15 and 43 DF, p-value: 1.076e-15
> 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.22522554 0.4504511 0.7747745
[2,] 0.28409259 0.5681852 0.7159074
[3,] 0.24432965 0.4886593 0.7556704
[4,] 0.30023292 0.6004658 0.6997671
[5,] 0.33451510 0.6690302 0.6654849
[6,] 0.31386891 0.6277378 0.6861311
[7,] 0.29392760 0.5878552 0.7060724
[8,] 0.26763920 0.5352784 0.7323608
[9,] 0.24565763 0.4913153 0.7543424
[10,] 0.21427239 0.4285448 0.7857276
[11,] 0.15348068 0.3069614 0.8465193
[12,] 0.11463940 0.2292788 0.8853606
[13,] 0.13978213 0.2795643 0.8602179
[14,] 0.10231067 0.2046213 0.8976893
[15,] 0.09869579 0.1973916 0.9013042
[16,] 0.05958247 0.1191649 0.9404175
[17,] 0.05160226 0.1032045 0.9483977
[18,] 0.15948163 0.3189633 0.8405184
[19,] 0.22162628 0.4432526 0.7783737
[20,] 0.14628414 0.2925683 0.8537159
[21,] 0.14220776 0.2844155 0.8577922
[22,] 0.08185896 0.1637179 0.9181410
> postscript(file="/var/www/html/rcomp/tmp/1xyt71258739139.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/2la461258739139.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/30dxe1258739139.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/4c0pv1258739139.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/5446g1258739139.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 = 59
Frequency = 1
1 2 3 4 5 6
-233.081488 -91.130399 -63.866899 -190.821524 94.691548 -107.941731
7 8 9 10 11 12
-126.703750 73.351233 75.401549 239.740498 126.933155 -19.428017
13 14 15 16 17 18
263.916354 72.965506 129.638313 -124.019656 -3.562356 -1.740132
19 20 21 22 23 24
-78.290032 -232.657155 -27.185760 -190.099173 -276.824180 -273.229078
25 26 27 28 29 30
-120.039824 28.866003 279.449466 217.251770 157.475372 229.003361
31 32 33 34 35 36
376.148098 -23.293381 -141.259790 187.441986 163.129265 508.717886
37 38 39 40 41 42
-121.950157 31.085390 57.680553 23.404062 -102.001954 -269.949123
43 44 45 46 47 48
-313.529486 176.579938 141.330791 -258.683462 240.999153 -216.060791
49 50 51 52 53 54
211.155114 -41.786501 -402.901433 74.185348 -146.602609 150.627625
55 56 57 58 59
142.375169 6.019364 -48.286790 21.600151 -254.237393
> postscript(file="/var/www/html/rcomp/tmp/6a8a91258739139.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 = 59
Frequency = 1
lag(myerror, k = 1) myerror
0 -233.081488 NA
1 -91.130399 -233.081488
2 -63.866899 -91.130399
3 -190.821524 -63.866899
4 94.691548 -190.821524
5 -107.941731 94.691548
6 -126.703750 -107.941731
7 73.351233 -126.703750
8 75.401549 73.351233
9 239.740498 75.401549
10 126.933155 239.740498
11 -19.428017 126.933155
12 263.916354 -19.428017
13 72.965506 263.916354
14 129.638313 72.965506
15 -124.019656 129.638313
16 -3.562356 -124.019656
17 -1.740132 -3.562356
18 -78.290032 -1.740132
19 -232.657155 -78.290032
20 -27.185760 -232.657155
21 -190.099173 -27.185760
22 -276.824180 -190.099173
23 -273.229078 -276.824180
24 -120.039824 -273.229078
25 28.866003 -120.039824
26 279.449466 28.866003
27 217.251770 279.449466
28 157.475372 217.251770
29 229.003361 157.475372
30 376.148098 229.003361
31 -23.293381 376.148098
32 -141.259790 -23.293381
33 187.441986 -141.259790
34 163.129265 187.441986
35 508.717886 163.129265
36 -121.950157 508.717886
37 31.085390 -121.950157
38 57.680553 31.085390
39 23.404062 57.680553
40 -102.001954 23.404062
41 -269.949123 -102.001954
42 -313.529486 -269.949123
43 176.579938 -313.529486
44 141.330791 176.579938
45 -258.683462 141.330791
46 240.999153 -258.683462
47 -216.060791 240.999153
48 211.155114 -216.060791
49 -41.786501 211.155114
50 -402.901433 -41.786501
51 74.185348 -402.901433
52 -146.602609 74.185348
53 150.627625 -146.602609
54 142.375169 150.627625
55 6.019364 142.375169
56 -48.286790 6.019364
57 21.600151 -48.286790
58 -254.237393 21.600151
59 NA -254.237393
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -91.130399 -233.081488
[2,] -63.866899 -91.130399
[3,] -190.821524 -63.866899
[4,] 94.691548 -190.821524
[5,] -107.941731 94.691548
[6,] -126.703750 -107.941731
[7,] 73.351233 -126.703750
[8,] 75.401549 73.351233
[9,] 239.740498 75.401549
[10,] 126.933155 239.740498
[11,] -19.428017 126.933155
[12,] 263.916354 -19.428017
[13,] 72.965506 263.916354
[14,] 129.638313 72.965506
[15,] -124.019656 129.638313
[16,] -3.562356 -124.019656
[17,] -1.740132 -3.562356
[18,] -78.290032 -1.740132
[19,] -232.657155 -78.290032
[20,] -27.185760 -232.657155
[21,] -190.099173 -27.185760
[22,] -276.824180 -190.099173
[23,] -273.229078 -276.824180
[24,] -120.039824 -273.229078
[25,] 28.866003 -120.039824
[26,] 279.449466 28.866003
[27,] 217.251770 279.449466
[28,] 157.475372 217.251770
[29,] 229.003361 157.475372
[30,] 376.148098 229.003361
[31,] -23.293381 376.148098
[32,] -141.259790 -23.293381
[33,] 187.441986 -141.259790
[34,] 163.129265 187.441986
[35,] 508.717886 163.129265
[36,] -121.950157 508.717886
[37,] 31.085390 -121.950157
[38,] 57.680553 31.085390
[39,] 23.404062 57.680553
[40,] -102.001954 23.404062
[41,] -269.949123 -102.001954
[42,] -313.529486 -269.949123
[43,] 176.579938 -313.529486
[44,] 141.330791 176.579938
[45,] -258.683462 141.330791
[46,] 240.999153 -258.683462
[47,] -216.060791 240.999153
[48,] 211.155114 -216.060791
[49,] -41.786501 211.155114
[50,] -402.901433 -41.786501
[51,] 74.185348 -402.901433
[52,] -146.602609 74.185348
[53,] 150.627625 -146.602609
[54,] 142.375169 150.627625
[55,] 6.019364 142.375169
[56,] -48.286790 6.019364
[57,] 21.600151 -48.286790
[58,] -254.237393 21.600151
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -91.130399 -233.081488
2 -63.866899 -91.130399
3 -190.821524 -63.866899
4 94.691548 -190.821524
5 -107.941731 94.691548
6 -126.703750 -107.941731
7 73.351233 -126.703750
8 75.401549 73.351233
9 239.740498 75.401549
10 126.933155 239.740498
11 -19.428017 126.933155
12 263.916354 -19.428017
13 72.965506 263.916354
14 129.638313 72.965506
15 -124.019656 129.638313
16 -3.562356 -124.019656
17 -1.740132 -3.562356
18 -78.290032 -1.740132
19 -232.657155 -78.290032
20 -27.185760 -232.657155
21 -190.099173 -27.185760
22 -276.824180 -190.099173
23 -273.229078 -276.824180
24 -120.039824 -273.229078
25 28.866003 -120.039824
26 279.449466 28.866003
27 217.251770 279.449466
28 157.475372 217.251770
29 229.003361 157.475372
30 376.148098 229.003361
31 -23.293381 376.148098
32 -141.259790 -23.293381
33 187.441986 -141.259790
34 163.129265 187.441986
35 508.717886 163.129265
36 -121.950157 508.717886
37 31.085390 -121.950157
38 57.680553 31.085390
39 23.404062 57.680553
40 -102.001954 23.404062
41 -269.949123 -102.001954
42 -313.529486 -269.949123
43 176.579938 -313.529486
44 141.330791 176.579938
45 -258.683462 141.330791
46 240.999153 -258.683462
47 -216.060791 240.999153
48 211.155114 -216.060791
49 -41.786501 211.155114
50 -402.901433 -41.786501
51 74.185348 -402.901433
52 -146.602609 74.185348
53 150.627625 -146.602609
54 142.375169 150.627625
55 6.019364 142.375169
56 -48.286790 6.019364
57 21.600151 -48.286790
58 -254.237393 21.600151
> 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/77jzx1258739139.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/87fi01258739139.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/9sfb91258739139.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/10pavu1258739139.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/11o4bn1258739139.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/12iuph1258739139.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/138rml1258739139.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/1422zb1258739139.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/15hk9f1258739139.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/16a4ui1258739139.tab")
+ }
>
> system("convert tmp/1xyt71258739139.ps tmp/1xyt71258739139.png")
> system("convert tmp/2la461258739139.ps tmp/2la461258739139.png")
> system("convert tmp/30dxe1258739139.ps tmp/30dxe1258739139.png")
> system("convert tmp/4c0pv1258739139.ps tmp/4c0pv1258739139.png")
> system("convert tmp/5446g1258739139.ps tmp/5446g1258739139.png")
> system("convert tmp/6a8a91258739139.ps tmp/6a8a91258739139.png")
> system("convert tmp/77jzx1258739139.ps tmp/77jzx1258739139.png")
> system("convert tmp/87fi01258739139.ps tmp/87fi01258739139.png")
> system("convert tmp/9sfb91258739139.ps tmp/9sfb91258739139.png")
> system("convert tmp/10pavu1258739139.ps tmp/10pavu1258739139.png")
>
>
> proc.time()
user system elapsed
2.452 1.588 3.471