R version 2.12.1 (2010-12-16)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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(9
+ ,20
+ ,1
+ ,14
+ ,3
+ ,1
+ ,1
+ ,9
+ ,14
+ ,1
+ ,8
+ ,3
+ ,0
+ ,1
+ ,9
+ ,18
+ ,0
+ ,12
+ ,6
+ ,1
+ ,1
+ ,9
+ ,12
+ ,1
+ ,7
+ ,2
+ ,0
+ ,1
+ ,9
+ ,16
+ ,0
+ ,10
+ ,1
+ ,1
+ ,0
+ ,9
+ ,13
+ ,0
+ ,7
+ ,2
+ ,0
+ ,0
+ ,9
+ ,22
+ ,1
+ ,16
+ ,8
+ ,1
+ ,1
+ ,9
+ ,16
+ ,1
+ ,11
+ ,1
+ ,1
+ ,0
+ ,9
+ ,20
+ ,0
+ ,14
+ ,4
+ ,1
+ ,1
+ ,9
+ ,10
+ ,0
+ ,6
+ ,0
+ ,0
+ ,0
+ ,9
+ ,22
+ ,0
+ ,16
+ ,4
+ ,1
+ ,0
+ ,9
+ ,17
+ ,1
+ ,11
+ ,2
+ ,0
+ ,1
+ ,9
+ ,21
+ ,0
+ ,16
+ ,1
+ ,1
+ ,1
+ ,9
+ ,18
+ ,1
+ ,12
+ ,2
+ ,1
+ ,1
+ ,9
+ ,13
+ ,0
+ ,7
+ ,3
+ ,0
+ ,0
+ ,9
+ ,17
+ ,0
+ ,13
+ ,1
+ ,1
+ ,0
+ ,9
+ ,17
+ ,1
+ ,11
+ ,2
+ ,1
+ ,1
+ ,9
+ ,19
+ ,1
+ ,15
+ ,6
+ ,1
+ ,0
+ ,9
+ ,12
+ ,1
+ ,7
+ ,0
+ ,0
+ ,1
+ ,9
+ ,14
+ ,1
+ ,9
+ ,1
+ ,0
+ ,1
+ ,9
+ ,13
+ ,0
+ ,7
+ ,3
+ ,0
+ ,1
+ ,9
+ ,20
+ ,1
+ ,14
+ ,5
+ ,1
+ ,1
+ ,9
+ ,20
+ ,1
+ ,15
+ ,0
+ ,1
+ ,1
+ ,9
+ ,13
+ ,1
+ ,7
+ ,1
+ ,0
+ ,1
+ ,9
+ ,21
+ ,1
+ ,15
+ ,3
+ ,1
+ ,1
+ ,9
+ ,21
+ ,1
+ ,17
+ ,6
+ ,1
+ ,1
+ ,9
+ ,19
+ ,1
+ ,15
+ ,5
+ ,1
+ ,0
+ ,9
+ ,18
+ ,1
+ ,14
+ ,4
+ ,1
+ ,0
+ ,9
+ ,20
+ ,0
+ ,14
+ ,4
+ ,0
+ ,0
+ ,9
+ ,14
+ ,1
+ ,8
+ ,4
+ ,1
+ ,1
+ ,9
+ ,14
+ ,0
+ ,8
+ ,0
+ ,0
+ ,1
+ ,9
+ ,20
+ ,1
+ ,14
+ ,3
+ ,1
+ ,0
+ ,9
+ ,21
+ ,1
+ ,14
+ ,5
+ ,1
+ ,1
+ ,9
+ ,14
+ ,0
+ ,8
+ ,3
+ ,0
+ ,0
+ ,9
+ ,16
+ ,1
+ ,11
+ ,1
+ ,1
+ ,1
+ ,9
+ ,21
+ ,1
+ ,16
+ ,5
+ ,1
+ ,1
+ ,9
+ ,16
+ ,1
+ ,10
+ ,5
+ ,1
+ ,1
+ ,9
+ ,14
+ ,1
+ ,8
+ ,0
+ ,0
+ ,1
+ ,9
+ ,19
+ ,1
+ ,14
+ ,3
+ ,1
+ ,1
+ ,9
+ ,22
+ ,1
+ ,16
+ ,6
+ ,1
+ ,0
+ ,9
+ ,19
+ ,0
+ ,13
+ ,3
+ ,1
+ ,1
+ ,9
+ ,11
+ ,1
+ ,5
+ ,1
+ ,0
+ ,0
+ ,9
+ ,13
+ ,1
+ ,8
+ ,2
+ ,0
+ ,1
+ ,9
+ ,16
+ ,1
+ ,10
+ ,2
+ ,0
+ ,0
+ ,9
+ ,14
+ ,0
+ ,8
+ ,2
+ ,0
+ ,1
+ ,9
+ ,19
+ ,1
+ ,13
+ ,4
+ ,1
+ ,1
+ ,9
+ ,21
+ ,1
+ ,15
+ ,4
+ ,1
+ ,1
+ ,9
+ ,12
+ ,0
+ ,6
+ ,0
+ ,0
+ ,1
+ ,9
+ ,17
+ ,0
+ ,12
+ ,3
+ ,1
+ ,1
+ ,9
+ ,21
+ ,1
+ ,16
+ ,6
+ ,0
+ ,1
+ ,9
+ ,11
+ ,1
+ ,5
+ ,3
+ ,1
+ ,0
+ ,9
+ ,19
+ ,0
+ ,15
+ ,1
+ ,1
+ ,1
+ ,9
+ ,18
+ ,0
+ ,12
+ ,4
+ ,1
+ ,0
+ ,9
+ ,14
+ ,0
+ ,8
+ ,3
+ ,0
+ ,1
+ ,9
+ ,19
+ ,0
+ ,13
+ ,3
+ ,1
+ ,1
+ ,9
+ ,20
+ ,1
+ ,14
+ ,3
+ ,1
+ ,1
+ ,10
+ ,18
+ ,0
+ ,12
+ ,2
+ ,1
+ ,1
+ ,10
+ ,22
+ ,0
+ ,16
+ ,6
+ ,1
+ ,1
+ ,10
+ ,16
+ ,1
+ ,10
+ ,5
+ ,1
+ ,1
+ ,10
+ ,20
+ ,0
+ ,15
+ ,5
+ ,1
+ ,0
+ ,10
+ ,14
+ ,0
+ ,8
+ ,2
+ ,0
+ ,1
+ ,10
+ ,22
+ ,1
+ ,16
+ ,4
+ ,1
+ ,1
+ ,10
+ ,25
+ ,0
+ ,19
+ ,2
+ ,1
+ ,1
+ ,10
+ ,20
+ ,0
+ ,14
+ ,5
+ ,1
+ ,0)
+ ,dim=c(7
+ ,64)
+ ,dimnames=list(c('Month'
+ ,'Income'
+ ,'Change'
+ ,'Size'
+ ,'Complex'
+ ,'Big4'
+ ,'Product')
+ ,1:64))
> y <- array(NA,dim=c(7,64),dimnames=list(c('Month','Income','Change','Size','Complex','Big4','Product'),1:64))
> 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 = 'Do not include Seasonal Dummies'
> par1 = '4'
> #'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
> 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
Size Month Income Change Complex Big4 Product t
1 14 9 20 1 3 1 1 1
2 8 9 14 1 3 0 1 2
3 12 9 18 0 6 1 1 3
4 7 9 12 1 2 0 1 4
5 10 9 16 0 1 1 0 5
6 7 9 13 0 2 0 0 6
7 16 9 22 1 8 1 1 7
8 11 9 16 1 1 1 0 8
9 14 9 20 0 4 1 1 9
10 6 9 10 0 0 0 0 10
11 16 9 22 0 4 1 0 11
12 11 9 17 1 2 0 1 12
13 16 9 21 0 1 1 1 13
14 12 9 18 1 2 1 1 14
15 7 9 13 0 3 0 0 15
16 13 9 17 0 1 1 0 16
17 11 9 17 1 2 1 1 17
18 15 9 19 1 6 1 0 18
19 7 9 12 1 0 0 1 19
20 9 9 14 1 1 0 1 20
21 7 9 13 0 3 0 1 21
22 14 9 20 1 5 1 1 22
23 15 9 20 1 0 1 1 23
24 7 9 13 1 1 0 1 24
25 15 9 21 1 3 1 1 25
26 17 9 21 1 6 1 1 26
27 15 9 19 1 5 1 0 27
28 14 9 18 1 4 1 0 28
29 14 9 20 0 4 0 0 29
30 8 9 14 1 4 1 1 30
31 8 9 14 0 0 0 1 31
32 14 9 20 1 3 1 0 32
33 14 9 21 1 5 1 1 33
34 8 9 14 0 3 0 0 34
35 11 9 16 1 1 1 1 35
36 16 9 21 1 5 1 1 36
37 10 9 16 1 5 1 1 37
38 8 9 14 1 0 0 1 38
39 14 9 19 1 3 1 1 39
40 16 9 22 1 6 1 0 40
41 13 9 19 0 3 1 1 41
42 5 9 11 1 1 0 0 42
43 8 9 13 1 2 0 1 43
44 10 9 16 1 2 0 0 44
45 8 9 14 0 2 0 1 45
46 13 9 19 1 4 1 1 46
47 15 9 21 1 4 1 1 47
48 6 9 12 0 0 0 1 48
49 12 9 17 0 3 1 1 49
50 16 9 21 1 6 0 1 50
51 5 9 11 1 3 1 0 51
52 15 9 19 0 1 1 1 52
53 12 9 18 0 4 1 0 53
54 8 9 14 0 3 0 1 54
55 13 9 19 0 3 1 1 55
56 14 9 20 1 3 1 1 56
57 12 10 18 0 2 1 1 57
58 16 10 22 0 6 1 1 58
59 10 10 16 1 5 1 1 59
60 15 10 20 0 5 1 0 60
61 8 10 14 0 2 0 1 61
62 16 10 22 1 4 1 1 62
63 19 10 25 0 2 1 1 63
64 14 10 20 0 5 1 0 64
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Month Income Change Complex Big4
-3.847735 -0.166870 1.000893 0.139650 -0.062965 0.253598
Product t
-0.306206 -0.002632
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-1.3545 -0.4191 -0.2534 0.4947 1.6900
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.847735 3.028356 -1.271 0.209
Month -0.166870 0.343715 -0.485 0.629
Income 1.000893 0.039990 25.029 <2e-16 ***
Change 0.139650 0.195672 0.714 0.478
Complex -0.062965 0.059605 -1.056 0.295
Big4 0.253598 0.262380 0.967 0.338
Product -0.306206 0.201329 -1.521 0.134
t -0.002632 0.005902 -0.446 0.657
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7111 on 56 degrees of freedom
Multiple R-squared: 0.9645, Adjusted R-squared: 0.96
F-statistic: 217.2 on 7 and 56 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.3208545 0.64170892 0.67914554
[2,] 0.3117561 0.62351220 0.68824390
[3,] 0.2723624 0.54472478 0.72763761
[4,] 0.4953013 0.99060263 0.50469869
[5,] 0.4789072 0.95781443 0.52109279
[6,] 0.5836771 0.83264589 0.41632295
[7,] 0.6209832 0.75803365 0.37901683
[8,] 0.7765791 0.44684177 0.22342088
[9,] 0.7075367 0.58492662 0.29246331
[10,] 0.6404250 0.71914992 0.35957496
[11,] 0.5979126 0.80417470 0.40208735
[12,] 0.5693769 0.86124620 0.43062310
[13,] 0.4862992 0.97259838 0.51370081
[14,] 0.4741344 0.94826871 0.52586564
[15,] 0.4460286 0.89205717 0.55397142
[16,] 0.6815106 0.63697886 0.31848943
[17,] 0.7545087 0.49098262 0.24549131
[18,] 0.8843292 0.23134170 0.11567085
[19,] 0.8540138 0.29197248 0.14598624
[20,] 0.8852804 0.22943917 0.11471958
[21,] 0.8460689 0.30786224 0.15393112
[22,] 0.8788854 0.24222918 0.12111459
[23,] 0.9639149 0.07217018 0.03608509
[24,] 0.9510940 0.09781210 0.04890605
[25,] 0.9378017 0.12439651 0.06219825
[26,] 0.9317139 0.13657222 0.06828611
[27,] 0.9094282 0.18114359 0.09057179
[28,] 0.8796769 0.24064622 0.12032311
[29,] 0.8706688 0.25866249 0.12933125
[30,] 0.8415480 0.31690402 0.15845201
[31,] 0.8133728 0.37325435 0.18662717
[32,] 0.7772903 0.44541936 0.22270968
[33,] 0.8050388 0.38992242 0.19496121
[34,] 0.7445645 0.51087092 0.25543546
[35,] 0.6869392 0.62612157 0.31306079
[36,] 0.6297179 0.74056424 0.37028212
[37,] 0.6535659 0.69286825 0.34643413
[38,] 0.6783608 0.64327849 0.32163924
[39,] 0.5893262 0.82134752 0.41067376
[40,] 0.5111701 0.97765985 0.48882992
[41,] 0.4585617 0.91712343 0.54143828
[42,] 0.9217013 0.15659744 0.07829872
[43,] 0.9422800 0.11544002 0.05772001
> postscript(file="/var/www/rcomp/tmp/1rwm41321900239.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/rcomp/tmp/2viol1321900239.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/rcomp/tmp/3881d1321900239.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/rcomp/tmp/4ay5n1321900239.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/rcomp/tmp/5lgk31321900239.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 = 64
Frequency = 1
1 2 3 4 5 6
-0.56381437 -0.30222467 -0.22821915 0.64186211 -0.84219711 -0.52032258
7 8 9 10 11 12
-0.23498322 0.02605022 -0.34014035 1.36695770 -0.64286826 -0.34154479
13 14 15 16 17 18
0.48060224 -0.59077087 -0.43366609 1.18586636 -0.58198031 1.36451783
19 20 21 22 23 24
0.55541932 0.61922983 -0.11166526 -0.38260410 0.30520523 -0.36934718
25 26 27 28 29 30
-0.50152932 1.68999698 1.32524508 1.26580616 -0.34010031 -0.41914967
31 32 33 34 35 36
-0.27512809 -0.78841528 -1.35454063 -0.38454318 0.40333208 0.65335666
37 38 39 40 41 42
-0.33954458 -0.39635105 0.53711124 -0.58024851 -0.31797394 -0.62638314
43 44 45 46 47 48
0.74363362 -0.56262000 -0.11234483 -0.38149713 -0.38065123 -0.22859025
49 50 51 52 53 54
0.70487204 1.00677298 -0.73035970 1.58505355 -0.52873313 -0.02568834
55 56 57 58 59 60
-0.28111991 -0.41903071 -0.17105627 0.07986158 -0.11476097 0.71774211
61 62 63 64
0.09664419 -0.17518791 -0.16151456 -0.27172817
> postscript(file="/var/www/rcomp/tmp/6q4u91321900239.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 = 64
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.56381437 NA
1 -0.30222467 -0.56381437
2 -0.22821915 -0.30222467
3 0.64186211 -0.22821915
4 -0.84219711 0.64186211
5 -0.52032258 -0.84219711
6 -0.23498322 -0.52032258
7 0.02605022 -0.23498322
8 -0.34014035 0.02605022
9 1.36695770 -0.34014035
10 -0.64286826 1.36695770
11 -0.34154479 -0.64286826
12 0.48060224 -0.34154479
13 -0.59077087 0.48060224
14 -0.43366609 -0.59077087
15 1.18586636 -0.43366609
16 -0.58198031 1.18586636
17 1.36451783 -0.58198031
18 0.55541932 1.36451783
19 0.61922983 0.55541932
20 -0.11166526 0.61922983
21 -0.38260410 -0.11166526
22 0.30520523 -0.38260410
23 -0.36934718 0.30520523
24 -0.50152932 -0.36934718
25 1.68999698 -0.50152932
26 1.32524508 1.68999698
27 1.26580616 1.32524508
28 -0.34010031 1.26580616
29 -0.41914967 -0.34010031
30 -0.27512809 -0.41914967
31 -0.78841528 -0.27512809
32 -1.35454063 -0.78841528
33 -0.38454318 -1.35454063
34 0.40333208 -0.38454318
35 0.65335666 0.40333208
36 -0.33954458 0.65335666
37 -0.39635105 -0.33954458
38 0.53711124 -0.39635105
39 -0.58024851 0.53711124
40 -0.31797394 -0.58024851
41 -0.62638314 -0.31797394
42 0.74363362 -0.62638314
43 -0.56262000 0.74363362
44 -0.11234483 -0.56262000
45 -0.38149713 -0.11234483
46 -0.38065123 -0.38149713
47 -0.22859025 -0.38065123
48 0.70487204 -0.22859025
49 1.00677298 0.70487204
50 -0.73035970 1.00677298
51 1.58505355 -0.73035970
52 -0.52873313 1.58505355
53 -0.02568834 -0.52873313
54 -0.28111991 -0.02568834
55 -0.41903071 -0.28111991
56 -0.17105627 -0.41903071
57 0.07986158 -0.17105627
58 -0.11476097 0.07986158
59 0.71774211 -0.11476097
60 0.09664419 0.71774211
61 -0.17518791 0.09664419
62 -0.16151456 -0.17518791
63 -0.27172817 -0.16151456
64 NA -0.27172817
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -0.30222467 -0.56381437
[2,] -0.22821915 -0.30222467
[3,] 0.64186211 -0.22821915
[4,] -0.84219711 0.64186211
[5,] -0.52032258 -0.84219711
[6,] -0.23498322 -0.52032258
[7,] 0.02605022 -0.23498322
[8,] -0.34014035 0.02605022
[9,] 1.36695770 -0.34014035
[10,] -0.64286826 1.36695770
[11,] -0.34154479 -0.64286826
[12,] 0.48060224 -0.34154479
[13,] -0.59077087 0.48060224
[14,] -0.43366609 -0.59077087
[15,] 1.18586636 -0.43366609
[16,] -0.58198031 1.18586636
[17,] 1.36451783 -0.58198031
[18,] 0.55541932 1.36451783
[19,] 0.61922983 0.55541932
[20,] -0.11166526 0.61922983
[21,] -0.38260410 -0.11166526
[22,] 0.30520523 -0.38260410
[23,] -0.36934718 0.30520523
[24,] -0.50152932 -0.36934718
[25,] 1.68999698 -0.50152932
[26,] 1.32524508 1.68999698
[27,] 1.26580616 1.32524508
[28,] -0.34010031 1.26580616
[29,] -0.41914967 -0.34010031
[30,] -0.27512809 -0.41914967
[31,] -0.78841528 -0.27512809
[32,] -1.35454063 -0.78841528
[33,] -0.38454318 -1.35454063
[34,] 0.40333208 -0.38454318
[35,] 0.65335666 0.40333208
[36,] -0.33954458 0.65335666
[37,] -0.39635105 -0.33954458
[38,] 0.53711124 -0.39635105
[39,] -0.58024851 0.53711124
[40,] -0.31797394 -0.58024851
[41,] -0.62638314 -0.31797394
[42,] 0.74363362 -0.62638314
[43,] -0.56262000 0.74363362
[44,] -0.11234483 -0.56262000
[45,] -0.38149713 -0.11234483
[46,] -0.38065123 -0.38149713
[47,] -0.22859025 -0.38065123
[48,] 0.70487204 -0.22859025
[49,] 1.00677298 0.70487204
[50,] -0.73035970 1.00677298
[51,] 1.58505355 -0.73035970
[52,] -0.52873313 1.58505355
[53,] -0.02568834 -0.52873313
[54,] -0.28111991 -0.02568834
[55,] -0.41903071 -0.28111991
[56,] -0.17105627 -0.41903071
[57,] 0.07986158 -0.17105627
[58,] -0.11476097 0.07986158
[59,] 0.71774211 -0.11476097
[60,] 0.09664419 0.71774211
[61,] -0.17518791 0.09664419
[62,] -0.16151456 -0.17518791
[63,] -0.27172817 -0.16151456
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -0.30222467 -0.56381437
2 -0.22821915 -0.30222467
3 0.64186211 -0.22821915
4 -0.84219711 0.64186211
5 -0.52032258 -0.84219711
6 -0.23498322 -0.52032258
7 0.02605022 -0.23498322
8 -0.34014035 0.02605022
9 1.36695770 -0.34014035
10 -0.64286826 1.36695770
11 -0.34154479 -0.64286826
12 0.48060224 -0.34154479
13 -0.59077087 0.48060224
14 -0.43366609 -0.59077087
15 1.18586636 -0.43366609
16 -0.58198031 1.18586636
17 1.36451783 -0.58198031
18 0.55541932 1.36451783
19 0.61922983 0.55541932
20 -0.11166526 0.61922983
21 -0.38260410 -0.11166526
22 0.30520523 -0.38260410
23 -0.36934718 0.30520523
24 -0.50152932 -0.36934718
25 1.68999698 -0.50152932
26 1.32524508 1.68999698
27 1.26580616 1.32524508
28 -0.34010031 1.26580616
29 -0.41914967 -0.34010031
30 -0.27512809 -0.41914967
31 -0.78841528 -0.27512809
32 -1.35454063 -0.78841528
33 -0.38454318 -1.35454063
34 0.40333208 -0.38454318
35 0.65335666 0.40333208
36 -0.33954458 0.65335666
37 -0.39635105 -0.33954458
38 0.53711124 -0.39635105
39 -0.58024851 0.53711124
40 -0.31797394 -0.58024851
41 -0.62638314 -0.31797394
42 0.74363362 -0.62638314
43 -0.56262000 0.74363362
44 -0.11234483 -0.56262000
45 -0.38149713 -0.11234483
46 -0.38065123 -0.38149713
47 -0.22859025 -0.38065123
48 0.70487204 -0.22859025
49 1.00677298 0.70487204
50 -0.73035970 1.00677298
51 1.58505355 -0.73035970
52 -0.52873313 1.58505355
53 -0.02568834 -0.52873313
54 -0.28111991 -0.02568834
55 -0.41903071 -0.28111991
56 -0.17105627 -0.41903071
57 0.07986158 -0.17105627
58 -0.11476097 0.07986158
59 0.71774211 -0.11476097
60 0.09664419 0.71774211
61 -0.17518791 0.09664419
62 -0.16151456 -0.17518791
63 -0.27172817 -0.16151456
> 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/rcomp/tmp/7nvjj1321900239.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/rcomp/tmp/8sldj1321900239.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/rcomp/tmp/9h0od1321900239.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/rcomp/tmp/10jczr1321900239.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/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/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/rcomp/tmp/11k15v1321900239.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/rcomp/tmp/120ajl1321900239.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/rcomp/tmp/13r3pi1321900239.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/rcomp/tmp/147jsv1321900239.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/rcomp/tmp/15fbe71321900239.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/rcomp/tmp/16wsre1321900239.tab")
+ }
>
> try(system("convert tmp/1rwm41321900239.ps tmp/1rwm41321900239.png",intern=TRUE))
character(0)
> try(system("convert tmp/2viol1321900239.ps tmp/2viol1321900239.png",intern=TRUE))
character(0)
> try(system("convert tmp/3881d1321900239.ps tmp/3881d1321900239.png",intern=TRUE))
character(0)
> try(system("convert tmp/4ay5n1321900239.ps tmp/4ay5n1321900239.png",intern=TRUE))
character(0)
> try(system("convert tmp/5lgk31321900239.ps tmp/5lgk31321900239.png",intern=TRUE))
character(0)
> try(system("convert tmp/6q4u91321900239.ps tmp/6q4u91321900239.png",intern=TRUE))
character(0)
> try(system("convert tmp/7nvjj1321900239.ps tmp/7nvjj1321900239.png",intern=TRUE))
character(0)
> try(system("convert tmp/8sldj1321900239.ps tmp/8sldj1321900239.png",intern=TRUE))
character(0)
> try(system("convert tmp/9h0od1321900239.ps tmp/9h0od1321900239.png",intern=TRUE))
character(0)
> try(system("convert tmp/10jczr1321900239.ps tmp/10jczr1321900239.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.132 0.740 4.849