R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i686-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(0
+ ,20
+ ,5
+ ,28
+ ,3
+ ,-2
+ ,23
+ ,6
+ ,24
+ ,1
+ ,-4
+ ,27
+ ,6
+ ,24
+ ,0
+ ,-6
+ ,23
+ ,6
+ ,28
+ ,1
+ ,-2
+ ,21
+ ,5
+ ,22
+ ,1
+ ,1
+ ,18
+ ,6
+ ,24
+ ,3
+ ,7
+ ,16
+ ,6
+ ,23
+ ,5
+ ,2
+ ,11
+ ,6
+ ,22
+ ,5
+ ,2
+ ,14
+ ,4
+ ,25
+ ,4
+ ,13
+ ,-3
+ ,6
+ ,23
+ ,11
+ ,7
+ ,2
+ ,5
+ ,21
+ ,8
+ ,-1
+ ,26
+ ,4
+ ,21
+ ,-1
+ ,1
+ ,11
+ ,6
+ ,19
+ ,4
+ ,0
+ ,11
+ ,4
+ ,21
+ ,4
+ ,0
+ ,11
+ ,6
+ ,23
+ ,4
+ ,5
+ ,3
+ ,6
+ ,16
+ ,6
+ ,3
+ ,8
+ ,5
+ ,22
+ ,6
+ ,6
+ ,8
+ ,5
+ ,20
+ ,6
+ ,7
+ ,7
+ ,4
+ ,19
+ ,6
+ ,-6
+ ,3
+ ,3
+ ,20
+ ,4
+ ,-8
+ ,4
+ ,2
+ ,14
+ ,1
+ ,-5
+ ,-7
+ ,4
+ ,19
+ ,6
+ ,-14
+ ,0
+ ,1
+ ,15
+ ,0
+ ,-13
+ ,-5
+ ,2
+ ,14
+ ,2
+ ,-15
+ ,5
+ ,-1
+ ,13
+ ,-2
+ ,-14
+ ,-1
+ ,2
+ ,11
+ ,0
+ ,-10
+ ,-4
+ ,0
+ ,11
+ ,1
+ ,-14
+ ,4
+ ,-1
+ ,9
+ ,-3
+ ,-18
+ ,7
+ ,0
+ ,12
+ ,-3
+ ,-22
+ ,6
+ ,-3
+ ,9
+ ,-5
+ ,-24
+ ,13
+ ,-2
+ ,11
+ ,-7
+ ,-17
+ ,20
+ ,-1
+ ,9
+ ,-7
+ ,-16
+ ,21
+ ,1
+ ,14
+ ,-5
+ ,-17
+ ,37
+ ,-5
+ ,8
+ ,-13
+ ,-22
+ ,52
+ ,-2
+ ,13
+ ,-16
+ ,-25
+ ,59
+ ,-4
+ ,8
+ ,-20
+ ,-18
+ ,66
+ ,-1
+ ,15
+ ,-18
+ ,-23
+ ,73
+ ,-1
+ ,12
+ ,-21
+ ,-20
+ ,71
+ ,-3
+ ,14
+ ,-20
+ ,-9
+ ,69
+ ,0
+ ,13
+ ,-16
+ ,-4
+ ,63
+ ,2
+ ,11
+ ,-14
+ ,0
+ ,68
+ ,2
+ ,16
+ ,-12
+ ,3
+ ,58
+ ,0
+ ,14
+ ,-10
+ ,14
+ ,50
+ ,3
+ ,19
+ ,-3
+ ,13
+ ,50
+ ,3
+ ,18
+ ,-4
+ ,12
+ ,50
+ ,4
+ ,16
+ ,-4
+ ,16
+ ,47
+ ,5
+ ,20
+ ,-1
+ ,7
+ ,60
+ ,3
+ ,17
+ ,-8
+ ,2
+ ,62
+ ,2
+ ,17
+ ,-10
+ ,1
+ ,63
+ ,1
+ ,18
+ ,-11
+ ,7
+ ,56
+ ,3
+ ,20
+ ,-7
+ ,10
+ ,38
+ ,3
+ ,17
+ ,-2
+ ,3
+ ,45
+ ,1
+ ,16
+ ,-6
+ ,2
+ ,39
+ ,3
+ ,16
+ ,-4
+ ,12
+ ,26
+ ,3
+ ,12
+ ,0
+ ,14
+ ,25
+ ,4
+ ,15
+ ,2
+ ,11
+ ,19
+ ,2
+ ,13
+ ,2
+ ,13
+ ,14
+ ,5
+ ,17
+ ,5
+ ,17
+ ,6
+ ,4
+ ,19
+ ,8
+ ,14
+ ,4
+ ,3
+ ,21
+ ,8
+ ,7
+ ,5
+ ,1
+ ,19
+ ,5
+ ,16
+ ,-3
+ ,4
+ ,20
+ ,10
+ ,5
+ ,-5
+ ,1
+ ,14
+ ,6
+ ,5
+ ,0
+ ,1
+ ,18
+ ,6
+ ,15
+ ,-6
+ ,3
+ ,14
+ ,9
+ ,9
+ ,4
+ ,1
+ ,15
+ ,5
+ ,4
+ ,-3
+ ,1
+ ,11
+ ,5
+ ,-9
+ ,14
+ ,2
+ ,6
+ ,-4
+ ,-14
+ ,16
+ ,0
+ ,11
+ ,-5
+ ,-4
+ ,17
+ ,3
+ ,13
+ ,-1
+ ,-19
+ ,25
+ ,0
+ ,14
+ ,-8
+ ,-10
+ ,25
+ ,-4
+ ,7
+ ,-8
+ ,-22
+ ,30
+ ,-2
+ ,1
+ ,-13
+ ,-25
+ ,51
+ ,-4
+ ,8
+ ,-18
+ ,-8
+ ,31
+ ,-1
+ ,8
+ ,-8
+ ,-8
+ ,31
+ ,-1
+ ,7
+ ,-8
+ ,-8
+ ,25
+ ,0
+ ,11
+ ,-6
+ ,-2
+ ,35
+ ,2
+ ,13
+ ,-5
+ ,-6
+ ,39
+ ,0
+ ,1
+ ,-11
+ ,-10
+ ,48
+ ,-1
+ ,4
+ ,-14
+ ,-11
+ ,41
+ ,0
+ ,4
+ ,-12
+ ,-14
+ ,47
+ ,-2
+ ,10
+ ,-13
+ ,-25
+ ,61
+ ,-1
+ ,8
+ ,-19)
+ ,dim=c(5
+ ,83)
+ ,dimnames=list(c('X_1'
+ ,'X_2'
+ ,'X_3'
+ ,'X_4'
+ ,'Y_1')
+ ,1:83))
> y <- array(NA,dim=c(5,83),dimnames=list(c('X_1','X_2','X_3','X_4','Y_1'),1:83))
> 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 = 'Do not include Seasonal Dummies'
> par1 = '1'
> par3 <- 'No Linear Trend'
> par2 <- 'Do not include Seasonal 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, 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
X_1 X_2 X_3 X_4 Y_1
1 0 20 5 28 3
2 -2 23 6 24 1
3 -4 27 6 24 0
4 -6 23 6 28 1
5 -2 21 5 22 1
6 1 18 6 24 3
7 7 16 6 23 5
8 2 11 6 22 5
9 2 14 4 25 4
10 13 -3 6 23 11
11 7 2 5 21 8
12 -1 26 4 21 -1
13 1 11 6 19 4
14 0 11 4 21 4
15 0 11 6 23 4
16 5 3 6 16 6
17 3 8 5 22 6
18 6 8 5 20 6
19 7 7 4 19 6
20 -6 3 3 20 4
21 -8 4 2 14 1
22 -5 -7 4 19 6
23 -14 0 1 15 0
24 -13 -5 2 14 2
25 -15 5 -1 13 -2
26 -14 -1 2 11 0
27 -10 -4 0 11 1
28 -14 4 -1 9 -3
29 -18 7 0 12 -3
30 -22 6 -3 9 -5
31 -24 13 -2 11 -7
32 -17 20 -1 9 -7
33 -16 21 1 14 -5
34 -17 37 -5 8 -13
35 -22 52 -2 13 -16
36 -25 59 -4 8 -20
37 -18 66 -1 15 -18
38 -23 73 -1 12 -21
39 -20 71 -3 14 -20
40 -9 69 0 13 -16
41 -4 63 2 11 -14
42 0 68 2 16 -12
43 3 58 0 14 -10
44 14 50 3 19 -3
45 13 50 3 18 -4
46 12 50 4 16 -4
47 16 47 5 20 -1
48 7 60 3 17 -8
49 2 62 2 17 -10
50 1 63 1 18 -11
51 7 56 3 20 -7
52 10 38 3 17 -2
53 3 45 1 16 -6
54 2 39 3 16 -4
55 12 26 3 12 0
56 14 25 4 15 2
57 11 19 2 13 2
58 13 14 5 17 5
59 17 6 4 19 8
60 14 4 3 21 8
61 7 5 1 19 5
62 16 -3 4 20 10
63 5 -5 1 14 6
64 5 0 1 18 6
65 15 -6 3 14 9
66 9 4 1 15 5
67 4 -3 1 11 5
68 -9 14 2 6 -4
69 -14 16 0 11 -5
70 -4 17 3 13 -1
71 -19 25 0 14 -8
72 -10 25 -4 7 -8
73 -22 30 -2 1 -13
74 -25 51 -4 8 -18
75 -8 31 -1 8 -8
76 -8 31 -1 7 -8
77 -8 25 0 11 -6
78 -2 35 2 13 -5
79 -6 39 0 1 -11
80 -10 48 -1 4 -14
81 -11 41 0 4 -12
82 -14 47 -2 10 -13
83 -25 61 -1 8 -19
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X_2 X_3 X_4 Y_1
-0.2937 0.9784 -1.1066 -0.9342 3.9768
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2.60386 -0.94119 0.02746 0.84619 2.82187
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.29368 0.57993 -0.506 0.614
X_2 0.97837 0.01868 52.388 < 2e-16 ***
X_3 -1.10657 0.11765 -9.406 1.77e-14 ***
X_4 -0.93425 0.04467 -20.914 < 2e-16 ***
Y_1 3.97683 0.07186 55.341 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.336 on 78 degrees of freedom
Multiple R-squared: 0.9879, Adjusted R-squared: 0.9872
F-statistic: 1586 on 4 and 78 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.3629936 0.7259873 0.6370064
[2,] 0.2097596 0.4195193 0.7902404
[3,] 0.1208693 0.2417387 0.8791307
[4,] 0.1079605 0.2159209 0.8920395
[5,] 0.2288121 0.4576241 0.7711879
[6,] 0.1507735 0.3015470 0.8492265
[7,] 0.3064611 0.6129222 0.6935389
[8,] 0.4283698 0.8567396 0.5716302
[9,] 0.3686807 0.7373614 0.6313193
[10,] 0.4680142 0.9360285 0.5319858
[11,] 0.3948753 0.7897506 0.6051247
[12,] 0.3273727 0.6547455 0.6726273
[13,] 0.3789225 0.7578450 0.6210775
[14,] 0.3843060 0.7686120 0.6156940
[15,] 0.3366030 0.6732059 0.6633970
[16,] 0.3435848 0.6871695 0.6564152
[17,] 0.2952516 0.5905033 0.7047484
[18,] 0.2373873 0.4747746 0.7626127
[19,] 0.1877974 0.3755949 0.8122026
[20,] 0.1734680 0.3469360 0.8265320
[21,] 0.2035032 0.4070063 0.7964968
[22,] 0.2487245 0.4974490 0.7512755
[23,] 0.3842861 0.7685722 0.6157139
[24,] 0.3384904 0.6769808 0.6615096
[25,] 0.3107722 0.6215444 0.6892278
[26,] 0.4753193 0.9506386 0.5246807
[27,] 0.5695811 0.8608378 0.4304189
[28,] 0.5156463 0.9687074 0.4843537
[29,] 0.4493492 0.8986985 0.5506508
[30,] 0.4917994 0.9835987 0.5082006
[31,] 0.4808357 0.9616713 0.5191643
[32,] 0.4217315 0.8434630 0.5782685
[33,] 0.3594519 0.7189038 0.6405481
[34,] 0.6073059 0.7853882 0.3926941
[35,] 0.5948614 0.8102772 0.4051386
[36,] 0.5303390 0.9393220 0.4696610
[37,] 0.5050876 0.9898249 0.4949124
[38,] 0.4793676 0.9587353 0.5206324
[39,] 0.4536904 0.9073809 0.5463096
[40,] 0.4701864 0.9403729 0.5298136
[41,] 0.4097700 0.8195400 0.5902300
[42,] 0.3556200 0.7112400 0.6443800
[43,] 0.3769794 0.7539588 0.6230206
[44,] 0.5905800 0.8188401 0.4094200
[45,] 0.5447110 0.9105780 0.4552890
[46,] 0.4762241 0.9524483 0.5237759
[47,] 0.5242095 0.9515809 0.4757905
[48,] 0.5668090 0.8663820 0.4331910
[49,] 0.5096581 0.9806839 0.4903419
[50,] 0.4843939 0.9687878 0.5156061
[51,] 0.4469101 0.8938202 0.5530899
[52,] 0.4819580 0.9639160 0.5180420
[53,] 0.4798593 0.9597186 0.5201407
[54,] 0.4650700 0.9301399 0.5349300
[55,] 0.6267456 0.7465089 0.3732544
[56,] 0.5495067 0.9009865 0.4504933
[57,] 0.4910781 0.9821562 0.5089219
[58,] 0.5172521 0.9654959 0.4827479
[59,] 0.4820481 0.9640963 0.5179519
[60,] 0.5137316 0.9725369 0.4862684
[61,] 0.4369393 0.8738786 0.5630607
[62,] 0.3454771 0.6909541 0.6545229
[63,] 0.3439279 0.6878559 0.6560721
[64,] 0.3471398 0.6942796 0.6528602
[65,] 0.2957558 0.5915117 0.7042442
[66,] 0.4706868 0.9413735 0.5293132
[67,] 0.3326436 0.6652872 0.6673564
[68,] 0.2159815 0.4319629 0.7840185
> postscript(file="/var/wessaorg/rcomp/tmp/1abxz1355341146.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/wessaorg/rcomp/tmp/20krg1355341146.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/wessaorg/rcomp/tmp/3p7h91355341146.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/wessaorg/rcomp/tmp/4gzun1355341146.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/wessaorg/rcomp/tmp/5lo4m1355341146.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 = 83
Frequency = 1
1 2 3 4 5 6
0.48762066 0.87574510 -1.06090504 0.61274206 -0.14258412 0.81393499
7 8 9 10 11 12
-0.11723477 -1.15963333 0.47169234 0.61081543 -1.32561355 1.87840574
13 14 15 16 17 18
-0.98555065 -2.33019421 1.75144631 0.08500193 -2.30792434 -1.17642281
19 20 21 22 23 24
-1.23887293 -2.54405337 -0.30399877 0.45830897 1.41399040 -0.47549794
25 26 27 28 29 30
-0.60584000 -0.23806541 0.50707257 1.61236358 -1.41342810 -2.60385793
31 32 33 34 35 36
-0.52371860 -1.13423701 -2.18187971 0.73391975 0.97981815 0.15416055
37 38 39 40 41 42
2.21136653 -0.50948095 0.12578537 -0.43933213 2.82187144 -1.35239384
43 44 45 46 47 48
-0.60399378 -1.62388623 0.41869493 -1.34323253 -1.49504534 -0.39193408
49 50 51 52 53 54
-0.50158457 1.32455391 2.34746378 0.27122653 -0.81743411 -1.68773205
55 56 57 58 59 60
1.38676117 0.32078925 -0.89063044 1.12743905 1.78583641 1.50450414
61 62 63 64 65 66
1.37498469 2.57175608 0.51060946 -0.64424427 1.77163043 0.61635787
67 68 69 70 71 72
-1.27204814 1.32245794 0.80065223 -0.89682797 1.72855992 -0.23746882
73 74 75 76 77 78
-0.63752092 0.02746085 0.14627265 -0.78797659 1.97215140 -1.70673984
79 80 81 82 83
0.81662909 1.63796576 0.63946694 -0.86157008 -2.45969707
> postscript(file="/var/wessaorg/rcomp/tmp/656io1355341146.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 = 83
Frequency = 1
lag(myerror, k = 1) myerror
0 0.48762066 NA
1 0.87574510 0.48762066
2 -1.06090504 0.87574510
3 0.61274206 -1.06090504
4 -0.14258412 0.61274206
5 0.81393499 -0.14258412
6 -0.11723477 0.81393499
7 -1.15963333 -0.11723477
8 0.47169234 -1.15963333
9 0.61081543 0.47169234
10 -1.32561355 0.61081543
11 1.87840574 -1.32561355
12 -0.98555065 1.87840574
13 -2.33019421 -0.98555065
14 1.75144631 -2.33019421
15 0.08500193 1.75144631
16 -2.30792434 0.08500193
17 -1.17642281 -2.30792434
18 -1.23887293 -1.17642281
19 -2.54405337 -1.23887293
20 -0.30399877 -2.54405337
21 0.45830897 -0.30399877
22 1.41399040 0.45830897
23 -0.47549794 1.41399040
24 -0.60584000 -0.47549794
25 -0.23806541 -0.60584000
26 0.50707257 -0.23806541
27 1.61236358 0.50707257
28 -1.41342810 1.61236358
29 -2.60385793 -1.41342810
30 -0.52371860 -2.60385793
31 -1.13423701 -0.52371860
32 -2.18187971 -1.13423701
33 0.73391975 -2.18187971
34 0.97981815 0.73391975
35 0.15416055 0.97981815
36 2.21136653 0.15416055
37 -0.50948095 2.21136653
38 0.12578537 -0.50948095
39 -0.43933213 0.12578537
40 2.82187144 -0.43933213
41 -1.35239384 2.82187144
42 -0.60399378 -1.35239384
43 -1.62388623 -0.60399378
44 0.41869493 -1.62388623
45 -1.34323253 0.41869493
46 -1.49504534 -1.34323253
47 -0.39193408 -1.49504534
48 -0.50158457 -0.39193408
49 1.32455391 -0.50158457
50 2.34746378 1.32455391
51 0.27122653 2.34746378
52 -0.81743411 0.27122653
53 -1.68773205 -0.81743411
54 1.38676117 -1.68773205
55 0.32078925 1.38676117
56 -0.89063044 0.32078925
57 1.12743905 -0.89063044
58 1.78583641 1.12743905
59 1.50450414 1.78583641
60 1.37498469 1.50450414
61 2.57175608 1.37498469
62 0.51060946 2.57175608
63 -0.64424427 0.51060946
64 1.77163043 -0.64424427
65 0.61635787 1.77163043
66 -1.27204814 0.61635787
67 1.32245794 -1.27204814
68 0.80065223 1.32245794
69 -0.89682797 0.80065223
70 1.72855992 -0.89682797
71 -0.23746882 1.72855992
72 -0.63752092 -0.23746882
73 0.02746085 -0.63752092
74 0.14627265 0.02746085
75 -0.78797659 0.14627265
76 1.97215140 -0.78797659
77 -1.70673984 1.97215140
78 0.81662909 -1.70673984
79 1.63796576 0.81662909
80 0.63946694 1.63796576
81 -0.86157008 0.63946694
82 -2.45969707 -0.86157008
83 NA -2.45969707
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 0.87574510 0.48762066
[2,] -1.06090504 0.87574510
[3,] 0.61274206 -1.06090504
[4,] -0.14258412 0.61274206
[5,] 0.81393499 -0.14258412
[6,] -0.11723477 0.81393499
[7,] -1.15963333 -0.11723477
[8,] 0.47169234 -1.15963333
[9,] 0.61081543 0.47169234
[10,] -1.32561355 0.61081543
[11,] 1.87840574 -1.32561355
[12,] -0.98555065 1.87840574
[13,] -2.33019421 -0.98555065
[14,] 1.75144631 -2.33019421
[15,] 0.08500193 1.75144631
[16,] -2.30792434 0.08500193
[17,] -1.17642281 -2.30792434
[18,] -1.23887293 -1.17642281
[19,] -2.54405337 -1.23887293
[20,] -0.30399877 -2.54405337
[21,] 0.45830897 -0.30399877
[22,] 1.41399040 0.45830897
[23,] -0.47549794 1.41399040
[24,] -0.60584000 -0.47549794
[25,] -0.23806541 -0.60584000
[26,] 0.50707257 -0.23806541
[27,] 1.61236358 0.50707257
[28,] -1.41342810 1.61236358
[29,] -2.60385793 -1.41342810
[30,] -0.52371860 -2.60385793
[31,] -1.13423701 -0.52371860
[32,] -2.18187971 -1.13423701
[33,] 0.73391975 -2.18187971
[34,] 0.97981815 0.73391975
[35,] 0.15416055 0.97981815
[36,] 2.21136653 0.15416055
[37,] -0.50948095 2.21136653
[38,] 0.12578537 -0.50948095
[39,] -0.43933213 0.12578537
[40,] 2.82187144 -0.43933213
[41,] -1.35239384 2.82187144
[42,] -0.60399378 -1.35239384
[43,] -1.62388623 -0.60399378
[44,] 0.41869493 -1.62388623
[45,] -1.34323253 0.41869493
[46,] -1.49504534 -1.34323253
[47,] -0.39193408 -1.49504534
[48,] -0.50158457 -0.39193408
[49,] 1.32455391 -0.50158457
[50,] 2.34746378 1.32455391
[51,] 0.27122653 2.34746378
[52,] -0.81743411 0.27122653
[53,] -1.68773205 -0.81743411
[54,] 1.38676117 -1.68773205
[55,] 0.32078925 1.38676117
[56,] -0.89063044 0.32078925
[57,] 1.12743905 -0.89063044
[58,] 1.78583641 1.12743905
[59,] 1.50450414 1.78583641
[60,] 1.37498469 1.50450414
[61,] 2.57175608 1.37498469
[62,] 0.51060946 2.57175608
[63,] -0.64424427 0.51060946
[64,] 1.77163043 -0.64424427
[65,] 0.61635787 1.77163043
[66,] -1.27204814 0.61635787
[67,] 1.32245794 -1.27204814
[68,] 0.80065223 1.32245794
[69,] -0.89682797 0.80065223
[70,] 1.72855992 -0.89682797
[71,] -0.23746882 1.72855992
[72,] -0.63752092 -0.23746882
[73,] 0.02746085 -0.63752092
[74,] 0.14627265 0.02746085
[75,] -0.78797659 0.14627265
[76,] 1.97215140 -0.78797659
[77,] -1.70673984 1.97215140
[78,] 0.81662909 -1.70673984
[79,] 1.63796576 0.81662909
[80,] 0.63946694 1.63796576
[81,] -0.86157008 0.63946694
[82,] -2.45969707 -0.86157008
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 0.87574510 0.48762066
2 -1.06090504 0.87574510
3 0.61274206 -1.06090504
4 -0.14258412 0.61274206
5 0.81393499 -0.14258412
6 -0.11723477 0.81393499
7 -1.15963333 -0.11723477
8 0.47169234 -1.15963333
9 0.61081543 0.47169234
10 -1.32561355 0.61081543
11 1.87840574 -1.32561355
12 -0.98555065 1.87840574
13 -2.33019421 -0.98555065
14 1.75144631 -2.33019421
15 0.08500193 1.75144631
16 -2.30792434 0.08500193
17 -1.17642281 -2.30792434
18 -1.23887293 -1.17642281
19 -2.54405337 -1.23887293
20 -0.30399877 -2.54405337
21 0.45830897 -0.30399877
22 1.41399040 0.45830897
23 -0.47549794 1.41399040
24 -0.60584000 -0.47549794
25 -0.23806541 -0.60584000
26 0.50707257 -0.23806541
27 1.61236358 0.50707257
28 -1.41342810 1.61236358
29 -2.60385793 -1.41342810
30 -0.52371860 -2.60385793
31 -1.13423701 -0.52371860
32 -2.18187971 -1.13423701
33 0.73391975 -2.18187971
34 0.97981815 0.73391975
35 0.15416055 0.97981815
36 2.21136653 0.15416055
37 -0.50948095 2.21136653
38 0.12578537 -0.50948095
39 -0.43933213 0.12578537
40 2.82187144 -0.43933213
41 -1.35239384 2.82187144
42 -0.60399378 -1.35239384
43 -1.62388623 -0.60399378
44 0.41869493 -1.62388623
45 -1.34323253 0.41869493
46 -1.49504534 -1.34323253
47 -0.39193408 -1.49504534
48 -0.50158457 -0.39193408
49 1.32455391 -0.50158457
50 2.34746378 1.32455391
51 0.27122653 2.34746378
52 -0.81743411 0.27122653
53 -1.68773205 -0.81743411
54 1.38676117 -1.68773205
55 0.32078925 1.38676117
56 -0.89063044 0.32078925
57 1.12743905 -0.89063044
58 1.78583641 1.12743905
59 1.50450414 1.78583641
60 1.37498469 1.50450414
61 2.57175608 1.37498469
62 0.51060946 2.57175608
63 -0.64424427 0.51060946
64 1.77163043 -0.64424427
65 0.61635787 1.77163043
66 -1.27204814 0.61635787
67 1.32245794 -1.27204814
68 0.80065223 1.32245794
69 -0.89682797 0.80065223
70 1.72855992 -0.89682797
71 -0.23746882 1.72855992
72 -0.63752092 -0.23746882
73 0.02746085 -0.63752092
74 0.14627265 0.02746085
75 -0.78797659 0.14627265
76 1.97215140 -0.78797659
77 -1.70673984 1.97215140
78 0.81662909 -1.70673984
79 1.63796576 0.81662909
80 0.63946694 1.63796576
81 -0.86157008 0.63946694
82 -2.45969707 -0.86157008
> 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/wessaorg/rcomp/tmp/74bh41355341146.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/wessaorg/rcomp/tmp/8zcz71355341146.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/wessaorg/rcomp/tmp/9fubu1355341146.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/wessaorg/rcomp/tmp/109cq61355341146.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11f7e01355341146.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/wessaorg/rcomp/tmp/12a3ix1355341146.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/wessaorg/rcomp/tmp/13vtw81355341146.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/wessaorg/rcomp/tmp/14jo991355341146.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/wessaorg/rcomp/tmp/15qlne1355341146.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/wessaorg/rcomp/tmp/16frnr1355341146.tab")
+ }
>
> try(system("convert tmp/1abxz1355341146.ps tmp/1abxz1355341146.png",intern=TRUE))
character(0)
> try(system("convert tmp/20krg1355341146.ps tmp/20krg1355341146.png",intern=TRUE))
character(0)
> try(system("convert tmp/3p7h91355341146.ps tmp/3p7h91355341146.png",intern=TRUE))
character(0)
> try(system("convert tmp/4gzun1355341146.ps tmp/4gzun1355341146.png",intern=TRUE))
character(0)
> try(system("convert tmp/5lo4m1355341146.ps tmp/5lo4m1355341146.png",intern=TRUE))
character(0)
> try(system("convert tmp/656io1355341146.ps tmp/656io1355341146.png",intern=TRUE))
character(0)
> try(system("convert tmp/74bh41355341146.ps tmp/74bh41355341146.png",intern=TRUE))
character(0)
> try(system("convert tmp/8zcz71355341146.ps tmp/8zcz71355341146.png",intern=TRUE))
character(0)
> try(system("convert tmp/9fubu1355341146.ps tmp/9fubu1355341146.png",intern=TRUE))
character(0)
> try(system("convert tmp/109cq61355341146.ps tmp/109cq61355341146.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
6.790 1.190 8.032