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(3613
+ ,12266.7
+ ,4286.1
+ ,3884.3
+ ,3142.7
+ ,3956.2
+ ,3730.5
+ ,12919.9
+ ,4348.1
+ ,3892.2
+ ,3884.3
+ ,3142.7
+ ,3481.3
+ ,11497.3
+ ,3949.3
+ ,3613
+ ,3892.2
+ ,3884.3
+ ,3649.5
+ ,12142
+ ,4166.7
+ ,3730.5
+ ,3613
+ ,3892.2
+ ,4215.2
+ ,13919.4
+ ,4217.9
+ ,3481.3
+ ,3730.5
+ ,3613
+ ,4066.6
+ ,12656.8
+ ,4528.2
+ ,3649.5
+ ,3481.3
+ ,3730.5
+ ,4196.8
+ ,12034.1
+ ,4232.2
+ ,4215.2
+ ,3649.5
+ ,3481.3
+ ,4536.6
+ ,13199.7
+ ,4470.9
+ ,4066.6
+ ,4215.2
+ ,3649.5
+ ,4441.6
+ ,10881.3
+ ,5121.2
+ ,4196.8
+ ,4066.6
+ ,4215.2
+ ,3548.3
+ ,11301.2
+ ,4170.8
+ ,4536.6
+ ,4196.8
+ ,4066.6
+ ,4735.9
+ ,13643.9
+ ,4398.6
+ ,4441.6
+ ,4536.6
+ ,4196.8
+ ,4130.6
+ ,12517
+ ,4491.4
+ ,3548.3
+ ,4441.6
+ ,4536.6
+ ,4356.2
+ ,13981.1
+ ,4251.8
+ ,4735.9
+ ,3548.3
+ ,4441.6
+ ,4159.6
+ ,14275.7
+ ,4901.9
+ ,4130.6
+ ,4735.9
+ ,3548.3
+ ,3988
+ ,13435
+ ,4745.2
+ ,4356.2
+ ,4130.6
+ ,4735.9
+ ,4167.8
+ ,13565.7
+ ,4666.9
+ ,4159.6
+ ,4356.2
+ ,4130.6
+ ,4902.2
+ ,16216.3
+ ,4210.4
+ ,3988
+ ,4159.6
+ ,4356.2
+ ,3909.4
+ ,12970
+ ,5273.6
+ ,4167.8
+ ,3988
+ ,4159.6
+ ,4697.6
+ ,14079.9
+ ,4095.3
+ ,4902.2
+ ,4167.8
+ ,3988
+ ,4308.9
+ ,14235
+ ,4610.1
+ ,3909.4
+ ,4902.2
+ ,4167.8
+ ,4420.4
+ ,12213.4
+ ,4718.1
+ ,4697.6
+ ,3909.4
+ ,4902.2
+ ,3544.2
+ ,12581
+ ,4185.5
+ ,4308.9
+ ,4697.6
+ ,3909.4
+ ,4433
+ ,14130.4
+ ,4314.7
+ ,4420.4
+ ,4308.9
+ ,4697.6
+ ,4479.7
+ ,14210.8
+ ,4422.6
+ ,3544.2
+ ,4420.4
+ ,4308.9
+ ,4533.2
+ ,14378.5
+ ,5059.2
+ ,4433
+ ,3544.2
+ ,4420.4
+ ,4237.5
+ ,13142.8
+ ,5043.6
+ ,4479.7
+ ,4433
+ ,3544.2
+ ,4207.4
+ ,13714.7
+ ,4436.6
+ ,4533.2
+ ,4479.7
+ ,4433
+ ,4394
+ ,13621.9
+ ,4922.6
+ ,4237.5
+ ,4533.2
+ ,4479.7
+ ,5148.4
+ ,15379.8
+ ,4454.8
+ ,4207.4
+ ,4237.5
+ ,4533.2
+ ,4202.2
+ ,13306.3
+ ,5058.7
+ ,4394
+ ,4207.4
+ ,4237.5
+ ,4682.5
+ ,14391.2
+ ,4768.9
+ ,5148.4
+ ,4394
+ ,4207.4
+ ,4884.3
+ ,14909.9
+ ,5171.8
+ ,4202.2
+ ,5148.4
+ ,4394
+ ,5288.9
+ ,14025.4
+ ,4989.3
+ ,4682.5
+ ,4202.2
+ ,5148.4
+ ,4505.2
+ ,12951.2
+ ,5202.1
+ ,4884.3
+ ,4682.5
+ ,4202.2
+ ,4611.5
+ ,14344.3
+ ,4838.4
+ ,5288.9
+ ,4884.3
+ ,4682.5
+ ,5104
+ ,16093.4
+ ,4876.5
+ ,4505.2
+ ,5288.9
+ ,4884.3
+ ,4586.6
+ ,15413.6
+ ,5875.5
+ ,4611.5
+ ,4505.2
+ ,5288.9
+ ,4529.3
+ ,14705.7
+ ,5717.9
+ ,5104
+ ,4611.5
+ ,4505.2
+ ,4504.1
+ ,15972.8
+ ,4778.8
+ ,4586.6
+ ,5104
+ ,4611.5
+ ,4604.9
+ ,16241.4
+ ,6195.9
+ ,4529.3
+ ,4586.6
+ ,5104
+ ,4795.4
+ ,16626.4
+ ,4625.4
+ ,4504.1
+ ,4529.3
+ ,4586.6
+ ,5391.1
+ ,17136.2
+ ,5549.8
+ ,4604.9
+ ,4504.1
+ ,4529.3
+ ,5213.9
+ ,15622.9
+ ,6397.6
+ ,4795.4
+ ,4604.9
+ ,4504.1
+ ,5415
+ ,18003.9
+ ,5856.7
+ ,5391.1
+ ,4795.4
+ ,4604.9
+ ,5990.3
+ ,16136.1
+ ,6343.8
+ ,5213.9
+ ,5391.1
+ ,4795.4
+ ,4241.8
+ ,14423.7
+ ,6615.5
+ ,5415
+ ,5213.9
+ ,5391.1
+ ,5677.6
+ ,16789.4
+ ,5904.6
+ ,5990.3
+ ,5415
+ ,5213.9
+ ,5164.2
+ ,16782.2
+ ,6861
+ ,4241.8
+ ,5990.3
+ ,5415
+ ,3962.3
+ ,14133.8
+ ,6553.5
+ ,5677.6
+ ,4241.8
+ ,5990.3
+ ,4011
+ ,12607
+ ,5481
+ ,5164.2
+ ,5677.6
+ ,4241.8
+ ,3310.3
+ ,12004.5
+ ,5435.3
+ ,3962.3
+ ,5164.2
+ ,5677.6
+ ,3837.3
+ ,12175.4
+ ,5278
+ ,4011
+ ,3962.3
+ ,5164.2
+ ,4145.3
+ ,13268
+ ,4671.8
+ ,3310.3
+ ,4011
+ ,3962.3
+ ,3796.7
+ ,12299.3
+ ,4891.5
+ ,3837.3
+ ,3310.3
+ ,4011
+ ,3849.6
+ ,11800.6
+ ,4241.6
+ ,4145.3
+ ,3837.3
+ ,3310.3
+ ,4285
+ ,13873.3
+ ,4152.1
+ ,3796.7
+ ,4145.3
+ ,3837.3
+ ,4189.6
+ ,12269.6
+ ,4484.4
+ ,3849.6
+ ,3796.7
+ ,4145.3)
+ ,dim=c(6
+ ,57)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Yt-1'
+ ,'Yt-2'
+ ,'Yt-3'
+ ,'Yt-4
')
+ ,1:57))
> y <- array(NA,dim=c(6,57),dimnames=list(c('Y','X','Yt-1','Yt-2','Yt-3','Yt-4
'),1:57))
> 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 Yt-1 Yt-2 Yt-3 Yt-4\r M1 M2 M3 M4 M5 M6 M7 M8 M9 M10
1 3613.0 12266.7 4286.1 3884.3 3142.7 3956.2 1 0 0 0 0 0 0 0 0 0
2 3730.5 12919.9 4348.1 3892.2 3884.3 3142.7 0 1 0 0 0 0 0 0 0 0
3 3481.3 11497.3 3949.3 3613.0 3892.2 3884.3 0 0 1 0 0 0 0 0 0 0
4 3649.5 12142.0 4166.7 3730.5 3613.0 3892.2 0 0 0 1 0 0 0 0 0 0
5 4215.2 13919.4 4217.9 3481.3 3730.5 3613.0 0 0 0 0 1 0 0 0 0 0
6 4066.6 12656.8 4528.2 3649.5 3481.3 3730.5 0 0 0 0 0 1 0 0 0 0
7 4196.8 12034.1 4232.2 4215.2 3649.5 3481.3 0 0 0 0 0 0 1 0 0 0
8 4536.6 13199.7 4470.9 4066.6 4215.2 3649.5 0 0 0 0 0 0 0 1 0 0
9 4441.6 10881.3 5121.2 4196.8 4066.6 4215.2 0 0 0 0 0 0 0 0 1 0
10 3548.3 11301.2 4170.8 4536.6 4196.8 4066.6 0 0 0 0 0 0 0 0 0 1
11 4735.9 13643.9 4398.6 4441.6 4536.6 4196.8 0 0 0 0 0 0 0 0 0 0
12 4130.6 12517.0 4491.4 3548.3 4441.6 4536.6 0 0 0 0 0 0 0 0 0 0
13 4356.2 13981.1 4251.8 4735.9 3548.3 4441.6 1 0 0 0 0 0 0 0 0 0
14 4159.6 14275.7 4901.9 4130.6 4735.9 3548.3 0 1 0 0 0 0 0 0 0 0
15 3988.0 13435.0 4745.2 4356.2 4130.6 4735.9 0 0 1 0 0 0 0 0 0 0
16 4167.8 13565.7 4666.9 4159.6 4356.2 4130.6 0 0 0 1 0 0 0 0 0 0
17 4902.2 16216.3 4210.4 3988.0 4159.6 4356.2 0 0 0 0 1 0 0 0 0 0
18 3909.4 12970.0 5273.6 4167.8 3988.0 4159.6 0 0 0 0 0 1 0 0 0 0
19 4697.6 14079.9 4095.3 4902.2 4167.8 3988.0 0 0 0 0 0 0 1 0 0 0
20 4308.9 14235.0 4610.1 3909.4 4902.2 4167.8 0 0 0 0 0 0 0 1 0 0
21 4420.4 12213.4 4718.1 4697.6 3909.4 4902.2 0 0 0 0 0 0 0 0 1 0
22 3544.2 12581.0 4185.5 4308.9 4697.6 3909.4 0 0 0 0 0 0 0 0 0 1
23 4433.0 14130.4 4314.7 4420.4 4308.9 4697.6 0 0 0 0 0 0 0 0 0 0
24 4479.7 14210.8 4422.6 3544.2 4420.4 4308.9 0 0 0 0 0 0 0 0 0 0
25 4533.2 14378.5 5059.2 4433.0 3544.2 4420.4 1 0 0 0 0 0 0 0 0 0
26 4237.5 13142.8 5043.6 4479.7 4433.0 3544.2 0 1 0 0 0 0 0 0 0 0
27 4207.4 13714.7 4436.6 4533.2 4479.7 4433.0 0 0 1 0 0 0 0 0 0 0
28 4394.0 13621.9 4922.6 4237.5 4533.2 4479.7 0 0 0 1 0 0 0 0 0 0
29 5148.4 15379.8 4454.8 4207.4 4237.5 4533.2 0 0 0 0 1 0 0 0 0 0
30 4202.2 13306.3 5058.7 4394.0 4207.4 4237.5 0 0 0 0 0 1 0 0 0 0
31 4682.5 14391.2 4768.9 5148.4 4394.0 4207.4 0 0 0 0 0 0 1 0 0 0
32 4884.3 14909.9 5171.8 4202.2 5148.4 4394.0 0 0 0 0 0 0 0 1 0 0
33 5288.9 14025.4 4989.3 4682.5 4202.2 5148.4 0 0 0 0 0 0 0 0 1 0
34 4505.2 12951.2 5202.1 4884.3 4682.5 4202.2 0 0 0 0 0 0 0 0 0 1
35 4611.5 14344.3 4838.4 5288.9 4884.3 4682.5 0 0 0 0 0 0 0 0 0 0
36 5104.0 16093.4 4876.5 4505.2 5288.9 4884.3 0 0 0 0 0 0 0 0 0 0
37 4586.6 15413.6 5875.5 4611.5 4505.2 5288.9 1 0 0 0 0 0 0 0 0 0
38 4529.3 14705.7 5717.9 5104.0 4611.5 4505.2 0 1 0 0 0 0 0 0 0 0
39 4504.1 15972.8 4778.8 4586.6 5104.0 4611.5 0 0 1 0 0 0 0 0 0 0
40 4604.9 16241.4 6195.9 4529.3 4586.6 5104.0 0 0 0 1 0 0 0 0 0 0
41 4795.4 16626.4 4625.4 4504.1 4529.3 4586.6 0 0 0 0 1 0 0 0 0 0
42 5391.1 17136.2 5549.8 4604.9 4504.1 4529.3 0 0 0 0 0 1 0 0 0 0
43 5213.9 15622.9 6397.6 4795.4 4604.9 4504.1 0 0 0 0 0 0 1 0 0 0
44 5415.0 18003.9 5856.7 5391.1 4795.4 4604.9 0 0 0 0 0 0 0 1 0 0
45 5990.3 16136.1 6343.8 5213.9 5391.1 4795.4 0 0 0 0 0 0 0 0 1 0
46 4241.8 14423.7 6615.5 5415.0 5213.9 5391.1 0 0 0 0 0 0 0 0 0 1
47 5677.6 16789.4 5904.6 5990.3 5415.0 5213.9 0 0 0 0 0 0 0 0 0 0
48 5164.2 16782.2 6861.0 4241.8 5990.3 5415.0 0 0 0 0 0 0 0 0 0 0
49 3962.3 14133.8 6553.5 5677.6 4241.8 5990.3 1 0 0 0 0 0 0 0 0 0
50 4011.0 12607.0 5481.0 5164.2 5677.6 4241.8 0 1 0 0 0 0 0 0 0 0
51 3310.3 12004.5 5435.3 3962.3 5164.2 5677.6 0 0 1 0 0 0 0 0 0 0
52 3837.3 12175.4 5278.0 4011.0 3962.3 5164.2 0 0 0 1 0 0 0 0 0 0
53 4145.3 13268.0 4671.8 3310.3 4011.0 3962.3 0 0 0 0 1 0 0 0 0 0
54 3796.7 12299.3 4891.5 3837.3 3310.3 4011.0 0 0 0 0 0 1 0 0 0 0
55 3849.6 11800.6 4241.6 4145.3 3837.3 3310.3 0 0 0 0 0 0 1 0 0 0
56 4285.0 13873.3 4152.1 3796.7 4145.3 3837.3 0 0 0 0 0 0 0 1 0 0
57 4189.6 12269.6 4484.4 3849.6 3796.7 4145.3 0 0 0 0 0 0 0 0 1 0
M11 t
1 0 1
2 0 2
3 0 3
4 0 4
5 0 5
6 0 6
7 0 7
8 0 8
9 0 9
10 0 10
11 1 11
12 0 12
13 0 13
14 0 14
15 0 15
16 0 16
17 0 17
18 0 18
19 0 19
20 0 20
21 0 21
22 0 22
23 1 23
24 0 24
25 0 25
26 0 26
27 0 27
28 0 28
29 0 29
30 0 30
31 0 31
32 0 32
33 0 33
34 0 34
35 1 35
36 0 36
37 0 37
38 0 38
39 0 39
40 0 40
41 0 41
42 0 42
43 0 43
44 0 44
45 0 45
46 0 46
47 1 47
48 0 48
49 0 49
50 0 50
51 0 51
52 0 52
53 0 53
54 0 54
55 0 55
56 0 56
57 0 57
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X `Yt-1` `Yt-2` `Yt-3` `Yt-4\r`
527.15396 0.24095 0.08186 0.21938 0.07324 -0.20323
M1 M2 M3 M4 M5 M6
-376.32012 -566.54873 -453.39293 -284.22635 -102.79622 -229.17430
M7 M8 M9 M10 M11 t
-101.47796 -162.03179 492.05456 -489.49492 -27.28533 -2.84688
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-329.93 -140.81 19.33 127.72 457.38
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 527.15396 428.00103 1.232 0.2254
X 0.24095 0.03207 7.513 4.28e-09 ***
`Yt-1` 0.08186 0.08442 0.970 0.3381
`Yt-2` 0.21938 0.11269 1.947 0.0588 .
`Yt-3` 0.07324 0.10921 0.671 0.5064
`Yt-4\r` -0.20323 0.12118 -1.677 0.1015
M1 -376.32012 223.34230 -1.685 0.1000 .
M2 -566.54873 229.55600 -2.468 0.0181 *
M3 -453.39293 163.56373 -2.772 0.0085 **
M4 -284.22635 175.82577 -1.617 0.1140
M5 -102.79622 172.23565 -0.597 0.5541
M6 -229.17430 207.46962 -1.105 0.2761
M7 -101.47796 234.64542 -0.432 0.6678
M8 -162.03179 175.80199 -0.922 0.3624
M9 492.05456 198.93774 2.473 0.0178 *
M10 -489.49492 222.37860 -2.201 0.0337 *
M11 -27.28533 202.41174 -0.135 0.8935
t -2.84688 2.32360 -1.225 0.2278
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 215.8 on 39 degrees of freedom
Multiple R-squared: 0.8984, Adjusted R-squared: 0.8542
F-statistic: 20.29 on 17 and 39 DF, p-value: 2.872e-14
> 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.10091573 0.20183145 0.8990843
[2,] 0.05431068 0.10862136 0.9456893
[3,] 0.07482784 0.14965567 0.9251722
[4,] 0.03379219 0.06758438 0.9662078
[5,] 0.11406806 0.22813613 0.8859319
[6,] 0.14576940 0.29153880 0.8542306
[7,] 0.08843856 0.17687713 0.9115614
[8,] 0.08905313 0.17810626 0.9109469
[9,] 0.16090955 0.32181910 0.8390904
[10,] 0.12771723 0.25543446 0.8722828
[11,] 0.11899831 0.23799662 0.8810017
[12,] 0.07978883 0.15957766 0.9202112
[13,] 0.08362731 0.16725462 0.9163727
[14,] 0.13765137 0.27530274 0.8623486
[15,] 0.29620484 0.59240968 0.7037952
[16,] 0.16995153 0.33990307 0.8300485
> postscript(file="/var/www/html/rcomp/tmp/1zmtz1258738769.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/2h9wl1258738769.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/3o3t51258738769.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/4a7qg1258738769.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/5e9zf1258738769.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 = 57
Frequency = 1
1 2 3 4 5 6
-119.862465 -193.128165 34.176209 -140.806058 -196.832223 67.852600
7 8 9 10 11 12
60.412349 188.565425 45.007844 -1.551177 165.966730 72.146651
13 14 15 16 17 18
129.334346 -134.129335 35.554574 -72.455575 -20.051900 -255.280325
19 20 21 22 23 24
27.939455 -176.325666 -188.740070 -299.740603 -90.023117 9.089147
25 26 27 28 29 30
241.087339 184.078523 121.028898 194.326953 423.997132 -41.619415
31 32 33 34 35 36
-109.138498 188.339688 287.013600 457.382129 -207.511254 19.280096
37 38 39 40 41 42
79.373783 123.525127 -141.384693 -237.076050 -284.788466 209.712669
43 44 45 46 47 48
148.602982 -240.480352 127.719560 -156.090348 131.567640 -100.515894
49 50 51 52 53 54
-329.933002 19.653851 -49.374987 256.010730 77.675457 19.334470
55 56 57
-127.816288 39.900905 -271.000933
> postscript(file="/var/www/html/rcomp/tmp/6j2lo1258738769.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 = 57
Frequency = 1
lag(myerror, k = 1) myerror
0 -119.862465 NA
1 -193.128165 -119.862465
2 34.176209 -193.128165
3 -140.806058 34.176209
4 -196.832223 -140.806058
5 67.852600 -196.832223
6 60.412349 67.852600
7 188.565425 60.412349
8 45.007844 188.565425
9 -1.551177 45.007844
10 165.966730 -1.551177
11 72.146651 165.966730
12 129.334346 72.146651
13 -134.129335 129.334346
14 35.554574 -134.129335
15 -72.455575 35.554574
16 -20.051900 -72.455575
17 -255.280325 -20.051900
18 27.939455 -255.280325
19 -176.325666 27.939455
20 -188.740070 -176.325666
21 -299.740603 -188.740070
22 -90.023117 -299.740603
23 9.089147 -90.023117
24 241.087339 9.089147
25 184.078523 241.087339
26 121.028898 184.078523
27 194.326953 121.028898
28 423.997132 194.326953
29 -41.619415 423.997132
30 -109.138498 -41.619415
31 188.339688 -109.138498
32 287.013600 188.339688
33 457.382129 287.013600
34 -207.511254 457.382129
35 19.280096 -207.511254
36 79.373783 19.280096
37 123.525127 79.373783
38 -141.384693 123.525127
39 -237.076050 -141.384693
40 -284.788466 -237.076050
41 209.712669 -284.788466
42 148.602982 209.712669
43 -240.480352 148.602982
44 127.719560 -240.480352
45 -156.090348 127.719560
46 131.567640 -156.090348
47 -100.515894 131.567640
48 -329.933002 -100.515894
49 19.653851 -329.933002
50 -49.374987 19.653851
51 256.010730 -49.374987
52 77.675457 256.010730
53 19.334470 77.675457
54 -127.816288 19.334470
55 39.900905 -127.816288
56 -271.000933 39.900905
57 NA -271.000933
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -193.128165 -119.862465
[2,] 34.176209 -193.128165
[3,] -140.806058 34.176209
[4,] -196.832223 -140.806058
[5,] 67.852600 -196.832223
[6,] 60.412349 67.852600
[7,] 188.565425 60.412349
[8,] 45.007844 188.565425
[9,] -1.551177 45.007844
[10,] 165.966730 -1.551177
[11,] 72.146651 165.966730
[12,] 129.334346 72.146651
[13,] -134.129335 129.334346
[14,] 35.554574 -134.129335
[15,] -72.455575 35.554574
[16,] -20.051900 -72.455575
[17,] -255.280325 -20.051900
[18,] 27.939455 -255.280325
[19,] -176.325666 27.939455
[20,] -188.740070 -176.325666
[21,] -299.740603 -188.740070
[22,] -90.023117 -299.740603
[23,] 9.089147 -90.023117
[24,] 241.087339 9.089147
[25,] 184.078523 241.087339
[26,] 121.028898 184.078523
[27,] 194.326953 121.028898
[28,] 423.997132 194.326953
[29,] -41.619415 423.997132
[30,] -109.138498 -41.619415
[31,] 188.339688 -109.138498
[32,] 287.013600 188.339688
[33,] 457.382129 287.013600
[34,] -207.511254 457.382129
[35,] 19.280096 -207.511254
[36,] 79.373783 19.280096
[37,] 123.525127 79.373783
[38,] -141.384693 123.525127
[39,] -237.076050 -141.384693
[40,] -284.788466 -237.076050
[41,] 209.712669 -284.788466
[42,] 148.602982 209.712669
[43,] -240.480352 148.602982
[44,] 127.719560 -240.480352
[45,] -156.090348 127.719560
[46,] 131.567640 -156.090348
[47,] -100.515894 131.567640
[48,] -329.933002 -100.515894
[49,] 19.653851 -329.933002
[50,] -49.374987 19.653851
[51,] 256.010730 -49.374987
[52,] 77.675457 256.010730
[53,] 19.334470 77.675457
[54,] -127.816288 19.334470
[55,] 39.900905 -127.816288
[56,] -271.000933 39.900905
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -193.128165 -119.862465
2 34.176209 -193.128165
3 -140.806058 34.176209
4 -196.832223 -140.806058
5 67.852600 -196.832223
6 60.412349 67.852600
7 188.565425 60.412349
8 45.007844 188.565425
9 -1.551177 45.007844
10 165.966730 -1.551177
11 72.146651 165.966730
12 129.334346 72.146651
13 -134.129335 129.334346
14 35.554574 -134.129335
15 -72.455575 35.554574
16 -20.051900 -72.455575
17 -255.280325 -20.051900
18 27.939455 -255.280325
19 -176.325666 27.939455
20 -188.740070 -176.325666
21 -299.740603 -188.740070
22 -90.023117 -299.740603
23 9.089147 -90.023117
24 241.087339 9.089147
25 184.078523 241.087339
26 121.028898 184.078523
27 194.326953 121.028898
28 423.997132 194.326953
29 -41.619415 423.997132
30 -109.138498 -41.619415
31 188.339688 -109.138498
32 287.013600 188.339688
33 457.382129 287.013600
34 -207.511254 457.382129
35 19.280096 -207.511254
36 79.373783 19.280096
37 123.525127 79.373783
38 -141.384693 123.525127
39 -237.076050 -141.384693
40 -284.788466 -237.076050
41 209.712669 -284.788466
42 148.602982 209.712669
43 -240.480352 148.602982
44 127.719560 -240.480352
45 -156.090348 127.719560
46 131.567640 -156.090348
47 -100.515894 131.567640
48 -329.933002 -100.515894
49 19.653851 -329.933002
50 -49.374987 19.653851
51 256.010730 -49.374987
52 77.675457 256.010730
53 19.334470 77.675457
54 -127.816288 19.334470
55 39.900905 -127.816288
56 -271.000933 39.900905
> 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/7fyit1258738769.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/88zc91258738769.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/9lolp1258738769.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/102m9a1258738769.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/11cm0c1258738769.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/12u17h1258738769.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/134rb81258738769.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/14z99n1258738769.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/15ybo61258738769.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/16lfe61258738769.tab")
+ }
>
> system("convert tmp/1zmtz1258738769.ps tmp/1zmtz1258738769.png")
> system("convert tmp/2h9wl1258738769.ps tmp/2h9wl1258738769.png")
> system("convert tmp/3o3t51258738769.ps tmp/3o3t51258738769.png")
> system("convert tmp/4a7qg1258738769.ps tmp/4a7qg1258738769.png")
> system("convert tmp/5e9zf1258738769.ps tmp/5e9zf1258738769.png")
> system("convert tmp/6j2lo1258738769.ps tmp/6j2lo1258738769.png")
> system("convert tmp/7fyit1258738769.ps tmp/7fyit1258738769.png")
> system("convert tmp/88zc91258738769.ps tmp/88zc91258738769.png")
> system("convert tmp/9lolp1258738769.ps tmp/9lolp1258738769.png")
> system("convert tmp/102m9a1258738769.ps tmp/102m9a1258738769.png")
>
>
> proc.time()
user system elapsed
2.388 1.576 2.756