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(149657,0,142773,0,133639,0,128332,0,120297,0,118632,0,155276,0,169316,0,167395,0,157939,0,149601,0,146310,0,141579,0,136473,0,129818,0,124226,0,116428,0,116440,0,147747,0,160069,0,163129,0,151108,0,141481,0,139174,0,134066,0,130104,0,123090,0,116598,0,109627,0,105428,0,137272,0,159836,0,155283,0,141514,0,131852,0,130691,0,128461,0,123066,0,117599,0,111599,0,105395,0,102334,0,131305,0,149033,0,144954,0,132404,0,122104,0,118755,0,116222,1,110924,1,103753,1,99983,1,93302,1,91496,1,119321,1,139261,1,133739,1,123913,1,113438,1,109416,1,109406,1,105645,1,101328,1,97686,1,93093,1,91382,1,122257,1,139183,1,139887,1,131822,1,116805,1,113706,1,113012,1,110452,1,107005,1,102841,1,98173,1,98181,1,137277,1,147579,1,146571,1,138920,1,130340,1,128140,1,127059,1,122860,1,117702,1,113537,1,108366,1,111078,1,150739,1,159129,1,157928,1,147768,1,137507,1,136919,1),dim=c(2,96),dimnames=list(c('Y','X'),1:96))
> y <- array(NA,dim=c(2,96),dimnames=list(c('Y','X'),1:96))
> for (i in 1:dim(x)[1])
+ {
+ for (j in 1:dim(x)[2])
+ {
+ y[i,j] <- as.numeric(x[i,j])
+ }
+ }
> par3 = 'No Linear Trend'
> par2 = 'Include Monthly Dummies'
> par1 = '1'
> #'GNU S' R Code compiled by R2WASP v. 1.0.44 ()
> #Author: Prof. Dr. P. Wessa
> #To cite this work: AUTHOR(S), (YEAR), YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/
> #Source of accompanying publication: Office for Research, Development, and Education
> #Technical description: Write here your technical program description (don't use hard returns!)
> library(lattice)
> library(lmtest)
Loading required package: zoo
Attaching package: 'zoo'
The following object(s) are masked from package:base :
as.Date.numeric
> n25 <- 25 #minimum number of obs. for Goldfeld-Quandt test
> par1 <- as.numeric(par1)
> x <- t(y)
> k <- length(x[1,])
> n <- length(x[,1])
> x1 <- cbind(x[,par1], x[,1:k!=par1])
> mycolnames <- c(colnames(x)[par1], colnames(x)[1:k!=par1])
> colnames(x1) <- mycolnames #colnames(x)[par1]
> x <- x1
> if (par3 == 'First Differences'){
+ x2 <- array(0, dim=c(n-1,k), dimnames=list(1:(n-1), paste('(1-B)',colnames(x),sep='')))
+ for (i in 1:n-1) {
+ for (j in 1:k) {
+ x2[i,j] <- x[i+1,j] - x[i,j]
+ }
+ }
+ x <- x2
+ }
> if (par2 == 'Include Monthly Dummies'){
+ x2 <- array(0, dim=c(n,11), dimnames=list(1:n, paste('M', seq(1:11), sep ='')))
+ for (i in 1:11){
+ x2[seq(i,n,12),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> if (par2 == 'Include Quarterly Dummies'){
+ x2 <- array(0, dim=c(n,3), dimnames=list(1:n, paste('Q', seq(1:3), sep ='')))
+ for (i in 1:3){
+ x2[seq(i,n,4),i] <- 1
+ }
+ x <- cbind(x, x2)
+ }
> k <- length(x[1,])
> if (par3 == 'Linear Trend'){
+ x <- cbind(x, c(1:n))
+ colnames(x)[k+1] <- 't'
+ }
> x
Y X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11
1 149657 0 1 0 0 0 0 0 0 0 0 0 0
2 142773 0 0 1 0 0 0 0 0 0 0 0 0
3 133639 0 0 0 1 0 0 0 0 0 0 0 0
4 128332 0 0 0 0 1 0 0 0 0 0 0 0
5 120297 0 0 0 0 0 1 0 0 0 0 0 0
6 118632 0 0 0 0 0 0 1 0 0 0 0 0
7 155276 0 0 0 0 0 0 0 1 0 0 0 0
8 169316 0 0 0 0 0 0 0 0 1 0 0 0
9 167395 0 0 0 0 0 0 0 0 0 1 0 0
10 157939 0 0 0 0 0 0 0 0 0 0 1 0
11 149601 0 0 0 0 0 0 0 0 0 0 0 1
12 146310 0 0 0 0 0 0 0 0 0 0 0 0
13 141579 0 1 0 0 0 0 0 0 0 0 0 0
14 136473 0 0 1 0 0 0 0 0 0 0 0 0
15 129818 0 0 0 1 0 0 0 0 0 0 0 0
16 124226 0 0 0 0 1 0 0 0 0 0 0 0
17 116428 0 0 0 0 0 1 0 0 0 0 0 0
18 116440 0 0 0 0 0 0 1 0 0 0 0 0
19 147747 0 0 0 0 0 0 0 1 0 0 0 0
20 160069 0 0 0 0 0 0 0 0 1 0 0 0
21 163129 0 0 0 0 0 0 0 0 0 1 0 0
22 151108 0 0 0 0 0 0 0 0 0 0 1 0
23 141481 0 0 0 0 0 0 0 0 0 0 0 1
24 139174 0 0 0 0 0 0 0 0 0 0 0 0
25 134066 0 1 0 0 0 0 0 0 0 0 0 0
26 130104 0 0 1 0 0 0 0 0 0 0 0 0
27 123090 0 0 0 1 0 0 0 0 0 0 0 0
28 116598 0 0 0 0 1 0 0 0 0 0 0 0
29 109627 0 0 0 0 0 1 0 0 0 0 0 0
30 105428 0 0 0 0 0 0 1 0 0 0 0 0
31 137272 0 0 0 0 0 0 0 1 0 0 0 0
32 159836 0 0 0 0 0 0 0 0 1 0 0 0
33 155283 0 0 0 0 0 0 0 0 0 1 0 0
34 141514 0 0 0 0 0 0 0 0 0 0 1 0
35 131852 0 0 0 0 0 0 0 0 0 0 0 1
36 130691 0 0 0 0 0 0 0 0 0 0 0 0
37 128461 0 1 0 0 0 0 0 0 0 0 0 0
38 123066 0 0 1 0 0 0 0 0 0 0 0 0
39 117599 0 0 0 1 0 0 0 0 0 0 0 0
40 111599 0 0 0 0 1 0 0 0 0 0 0 0
41 105395 0 0 0 0 0 1 0 0 0 0 0 0
42 102334 0 0 0 0 0 0 1 0 0 0 0 0
43 131305 0 0 0 0 0 0 0 1 0 0 0 0
44 149033 0 0 0 0 0 0 0 0 1 0 0 0
45 144954 0 0 0 0 0 0 0 0 0 1 0 0
46 132404 0 0 0 0 0 0 0 0 0 0 1 0
47 122104 0 0 0 0 0 0 0 0 0 0 0 1
48 118755 0 0 0 0 0 0 0 0 0 0 0 0
49 116222 1 1 0 0 0 0 0 0 0 0 0 0
50 110924 1 0 1 0 0 0 0 0 0 0 0 0
51 103753 1 0 0 1 0 0 0 0 0 0 0 0
52 99983 1 0 0 0 1 0 0 0 0 0 0 0
53 93302 1 0 0 0 0 1 0 0 0 0 0 0
54 91496 1 0 0 0 0 0 1 0 0 0 0 0
55 119321 1 0 0 0 0 0 0 1 0 0 0 0
56 139261 1 0 0 0 0 0 0 0 1 0 0 0
57 133739 1 0 0 0 0 0 0 0 0 1 0 0
58 123913 1 0 0 0 0 0 0 0 0 0 1 0
59 113438 1 0 0 0 0 0 0 0 0 0 0 1
60 109416 1 0 0 0 0 0 0 0 0 0 0 0
61 109406 1 1 0 0 0 0 0 0 0 0 0 0
62 105645 1 0 1 0 0 0 0 0 0 0 0 0
63 101328 1 0 0 1 0 0 0 0 0 0 0 0
64 97686 1 0 0 0 1 0 0 0 0 0 0 0
65 93093 1 0 0 0 0 1 0 0 0 0 0 0
66 91382 1 0 0 0 0 0 1 0 0 0 0 0
67 122257 1 0 0 0 0 0 0 1 0 0 0 0
68 139183 1 0 0 0 0 0 0 0 1 0 0 0
69 139887 1 0 0 0 0 0 0 0 0 1 0 0
70 131822 1 0 0 0 0 0 0 0 0 0 1 0
71 116805 1 0 0 0 0 0 0 0 0 0 0 1
72 113706 1 0 0 0 0 0 0 0 0 0 0 0
73 113012 1 1 0 0 0 0 0 0 0 0 0 0
74 110452 1 0 1 0 0 0 0 0 0 0 0 0
75 107005 1 0 0 1 0 0 0 0 0 0 0 0
76 102841 1 0 0 0 1 0 0 0 0 0 0 0
77 98173 1 0 0 0 0 1 0 0 0 0 0 0
78 98181 1 0 0 0 0 0 1 0 0 0 0 0
79 137277 1 0 0 0 0 0 0 1 0 0 0 0
80 147579 1 0 0 0 0 0 0 0 1 0 0 0
81 146571 1 0 0 0 0 0 0 0 0 1 0 0
82 138920 1 0 0 0 0 0 0 0 0 0 1 0
83 130340 1 0 0 0 0 0 0 0 0 0 0 1
84 128140 1 0 0 0 0 0 0 0 0 0 0 0
85 127059 1 1 0 0 0 0 0 0 0 0 0 0
86 122860 1 0 1 0 0 0 0 0 0 0 0 0
87 117702 1 0 0 1 0 0 0 0 0 0 0 0
88 113537 1 0 0 0 1 0 0 0 0 0 0 0
89 108366 1 0 0 0 0 1 0 0 0 0 0 0
90 111078 1 0 0 0 0 0 1 0 0 0 0 0
91 150739 1 0 0 0 0 0 0 1 0 0 0 0
92 159129 1 0 0 0 0 0 0 0 1 0 0 0
93 157928 1 0 0 0 0 0 0 0 0 1 0 0
94 147768 1 0 0 0 0 0 0 0 0 0 1 0
95 137507 1 0 0 0 0 0 0 0 0 0 0 1
96 136919 1 0 0 0 0 0 0 0 0 0 0 0
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X M1 M2 M3 M4
135213.4 -14649.1 -456.1 -5101.7 -11147.1 -16038.6
M5 M6 M7 M8 M9 M10
-22303.8 -23517.5 9760.4 25036.9 23221.9 12784.6
M11
2502.1
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-16458 -6472 -1252 7033 20414
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 135213.4 3374.1 40.075 < 2e-16 ***
X -14649.1 1871.6 -7.827 1.44e-11 ***
M1 -456.1 4584.4 -0.099 0.920986
M2 -5101.7 4584.4 -1.113 0.268989
M3 -11147.1 4584.4 -2.432 0.017188 *
M4 -16038.6 4584.4 -3.498 0.000755 ***
M5 -22303.8 4584.4 -4.865 5.37e-06 ***
M6 -23517.5 4584.4 -5.130 1.87e-06 ***
M7 9760.4 4584.4 2.129 0.036216 *
M8 25036.9 4584.4 5.461 4.83e-07 ***
M9 23221.9 4584.4 5.065 2.42e-06 ***
M10 12784.6 4584.4 2.789 0.006561 **
M11 2502.1 4584.4 0.546 0.586676
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9169 on 83 degrees of freedom
Multiple R-squared: 0.8019, Adjusted R-squared: 0.7733
F-statistic: 28 on 12 and 83 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.18732375 0.37464750 0.8126763
[2,] 0.09981834 0.19963668 0.9001817
[3,] 0.04660669 0.09321338 0.9533933
[4,] 0.04146907 0.08293814 0.9585309
[5,] 0.04629994 0.09259988 0.9537001
[6,] 0.02822395 0.05644790 0.9717761
[7,] 0.02240747 0.04481494 0.9775925
[8,] 0.02165853 0.04331706 0.9783415
[9,] 0.01855062 0.03710125 0.9814494
[10,] 0.03716908 0.07433817 0.9628309
[11,] 0.04570710 0.09141420 0.9542929
[12,] 0.04775996 0.09551992 0.9522400
[13,] 0.05353299 0.10706599 0.9464670
[14,] 0.05257665 0.10515330 0.9474233
[15,] 0.06930435 0.13860871 0.9306956
[16,] 0.10369840 0.20739680 0.8963016
[17,] 0.08381571 0.16763143 0.9161843
[18,] 0.08595584 0.17191167 0.9140442
[19,] 0.10421874 0.20843747 0.8957813
[20,] 0.12773516 0.25547033 0.8722648
[21,] 0.13830416 0.27660831 0.8616958
[22,] 0.16967984 0.33935969 0.8303202
[23,] 0.19948622 0.39897245 0.8005138
[24,] 0.20858867 0.41717734 0.7914113
[25,] 0.21470468 0.42940936 0.7852953
[26,] 0.20857809 0.41715618 0.7914219
[27,] 0.20516738 0.41033476 0.7948326
[28,] 0.23194521 0.46389042 0.7680548
[29,] 0.24510575 0.49021149 0.7548943
[30,] 0.28065823 0.56131645 0.7193418
[31,] 0.31709195 0.63418390 0.6829080
[32,] 0.35854802 0.71709604 0.6414520
[33,] 0.40374398 0.80748795 0.5962560
[34,] 0.33918017 0.67836033 0.6608198
[35,] 0.27990991 0.55981982 0.7200901
[36,] 0.23000864 0.46001728 0.7699914
[37,] 0.18572755 0.37145510 0.8142725
[38,] 0.15035727 0.30071455 0.8496427
[39,] 0.12415116 0.24830232 0.8758488
[40,] 0.13246042 0.26492085 0.8675396
[41,] 0.11034131 0.22068262 0.8896587
[42,] 0.10736409 0.21472818 0.8926359
[43,] 0.11016415 0.22032829 0.8898359
[44,] 0.11128252 0.22256505 0.8887175
[45,] 0.12691844 0.25383687 0.8730816
[46,] 0.11671569 0.23343139 0.8832843
[47,] 0.10549896 0.21099791 0.8945010
[48,] 0.09280385 0.18560770 0.9071961
[49,] 0.07946785 0.15893569 0.9205322
[50,] 0.06735047 0.13470094 0.9326495
[51,] 0.06541403 0.13082807 0.9345860
[52,] 0.11666753 0.23333506 0.8833325
[53,] 0.12588130 0.25176260 0.8741187
[54,] 0.12761075 0.25522150 0.8723892
[55,] 0.13008949 0.26017898 0.8699105
[56,] 0.18400166 0.36800331 0.8159983
[57,] 0.30998303 0.61996606 0.6900170
[58,] 0.33532221 0.67064442 0.6646778
[59,] 0.33933978 0.67867956 0.6606602
[60,] 0.32161511 0.64323022 0.6783849
[61,] 0.30493931 0.60987861 0.6950607
[62,] 0.28293932 0.56587864 0.7170607
[63,] 0.31215431 0.62430863 0.6878457
[64,] 0.38847775 0.77695550 0.6115222
[65,] 0.39815401 0.79630803 0.6018460
> postscript(file="/var/www/html/rcomp/tmp/1m7291258545954.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/2s6a41258545954.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/3gl0n1258545954.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/45rm91258545954.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/5rxzx1258545955.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 = 96
Frequency = 1
1 2 3 4 5 6
14899.6875 12661.3125 9572.6875 9157.1875 7387.3125 6936.0625
7 8 9 10 11 12
10302.1875 9065.6875 8959.6875 9940.9375 11885.4375 11096.5625
13 14 15 16 17 18
6821.6875 6361.3125 5751.6875 5051.1875 3518.3125 4744.0625
19 20 21 22 23 24
2773.1875 -181.3125 4693.6875 3109.9375 3765.4375 3960.5625
25 26 27 28 29 30
-691.3125 -7.6875 -976.3125 -2576.8125 -3282.6875 -6267.9375
31 32 33 34 35 36
-7701.8125 -414.3125 -3152.3125 -6484.0625 -5863.5625 -4522.4375
37 38 39 40 41 42
-6296.3125 -7045.6875 -6467.3125 -7575.8125 -7514.6875 -9361.9375
43 44 45 46 47 48
-13668.8125 -11217.3125 -13481.3125 -15594.0625 -15611.5625 -16458.4375
49 50 51 52 53 54
-3886.1875 -4538.5625 -5664.1875 -4542.6875 -4958.5625 -5550.8125
55 56 57 58 59 60
-11003.6875 -6340.1875 -10047.1875 -9435.9375 -9628.4375 -11148.3125
61 62 63 64 65 66
-10702.1875 -9817.5625 -8089.1875 -6839.6875 -5167.5625 -5664.8125
67 68 69 70 71 72
-8067.6875 -6418.1875 -3899.1875 -1526.9375 -6261.4375 -6858.3125
73 74 75 76 77 78
-7096.1875 -5010.5625 -2412.1875 -1684.6875 -87.5625 1134.1875
79 80 81 82 83 84
6952.3125 1977.8125 2784.8125 5571.0625 7273.5625 7575.6875
85 86 87 88 89 90
6950.8125 7397.4375 8284.8125 9011.3125 10105.4375 14031.1875
91 92 93 94 95 96
20414.3125 13527.8125 14141.8125 14419.0625 14440.5625 16354.6875
> postscript(file="/var/www/html/rcomp/tmp/60vup1258545955.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 = 96
Frequency = 1
lag(myerror, k = 1) myerror
0 14899.6875 NA
1 12661.3125 14899.6875
2 9572.6875 12661.3125
3 9157.1875 9572.6875
4 7387.3125 9157.1875
5 6936.0625 7387.3125
6 10302.1875 6936.0625
7 9065.6875 10302.1875
8 8959.6875 9065.6875
9 9940.9375 8959.6875
10 11885.4375 9940.9375
11 11096.5625 11885.4375
12 6821.6875 11096.5625
13 6361.3125 6821.6875
14 5751.6875 6361.3125
15 5051.1875 5751.6875
16 3518.3125 5051.1875
17 4744.0625 3518.3125
18 2773.1875 4744.0625
19 -181.3125 2773.1875
20 4693.6875 -181.3125
21 3109.9375 4693.6875
22 3765.4375 3109.9375
23 3960.5625 3765.4375
24 -691.3125 3960.5625
25 -7.6875 -691.3125
26 -976.3125 -7.6875
27 -2576.8125 -976.3125
28 -3282.6875 -2576.8125
29 -6267.9375 -3282.6875
30 -7701.8125 -6267.9375
31 -414.3125 -7701.8125
32 -3152.3125 -414.3125
33 -6484.0625 -3152.3125
34 -5863.5625 -6484.0625
35 -4522.4375 -5863.5625
36 -6296.3125 -4522.4375
37 -7045.6875 -6296.3125
38 -6467.3125 -7045.6875
39 -7575.8125 -6467.3125
40 -7514.6875 -7575.8125
41 -9361.9375 -7514.6875
42 -13668.8125 -9361.9375
43 -11217.3125 -13668.8125
44 -13481.3125 -11217.3125
45 -15594.0625 -13481.3125
46 -15611.5625 -15594.0625
47 -16458.4375 -15611.5625
48 -3886.1875 -16458.4375
49 -4538.5625 -3886.1875
50 -5664.1875 -4538.5625
51 -4542.6875 -5664.1875
52 -4958.5625 -4542.6875
53 -5550.8125 -4958.5625
54 -11003.6875 -5550.8125
55 -6340.1875 -11003.6875
56 -10047.1875 -6340.1875
57 -9435.9375 -10047.1875
58 -9628.4375 -9435.9375
59 -11148.3125 -9628.4375
60 -10702.1875 -11148.3125
61 -9817.5625 -10702.1875
62 -8089.1875 -9817.5625
63 -6839.6875 -8089.1875
64 -5167.5625 -6839.6875
65 -5664.8125 -5167.5625
66 -8067.6875 -5664.8125
67 -6418.1875 -8067.6875
68 -3899.1875 -6418.1875
69 -1526.9375 -3899.1875
70 -6261.4375 -1526.9375
71 -6858.3125 -6261.4375
72 -7096.1875 -6858.3125
73 -5010.5625 -7096.1875
74 -2412.1875 -5010.5625
75 -1684.6875 -2412.1875
76 -87.5625 -1684.6875
77 1134.1875 -87.5625
78 6952.3125 1134.1875
79 1977.8125 6952.3125
80 2784.8125 1977.8125
81 5571.0625 2784.8125
82 7273.5625 5571.0625
83 7575.6875 7273.5625
84 6950.8125 7575.6875
85 7397.4375 6950.8125
86 8284.8125 7397.4375
87 9011.3125 8284.8125
88 10105.4375 9011.3125
89 14031.1875 10105.4375
90 20414.3125 14031.1875
91 13527.8125 20414.3125
92 14141.8125 13527.8125
93 14419.0625 14141.8125
94 14440.5625 14419.0625
95 16354.6875 14440.5625
96 NA 16354.6875
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 12661.3125 14899.6875
[2,] 9572.6875 12661.3125
[3,] 9157.1875 9572.6875
[4,] 7387.3125 9157.1875
[5,] 6936.0625 7387.3125
[6,] 10302.1875 6936.0625
[7,] 9065.6875 10302.1875
[8,] 8959.6875 9065.6875
[9,] 9940.9375 8959.6875
[10,] 11885.4375 9940.9375
[11,] 11096.5625 11885.4375
[12,] 6821.6875 11096.5625
[13,] 6361.3125 6821.6875
[14,] 5751.6875 6361.3125
[15,] 5051.1875 5751.6875
[16,] 3518.3125 5051.1875
[17,] 4744.0625 3518.3125
[18,] 2773.1875 4744.0625
[19,] -181.3125 2773.1875
[20,] 4693.6875 -181.3125
[21,] 3109.9375 4693.6875
[22,] 3765.4375 3109.9375
[23,] 3960.5625 3765.4375
[24,] -691.3125 3960.5625
[25,] -7.6875 -691.3125
[26,] -976.3125 -7.6875
[27,] -2576.8125 -976.3125
[28,] -3282.6875 -2576.8125
[29,] -6267.9375 -3282.6875
[30,] -7701.8125 -6267.9375
[31,] -414.3125 -7701.8125
[32,] -3152.3125 -414.3125
[33,] -6484.0625 -3152.3125
[34,] -5863.5625 -6484.0625
[35,] -4522.4375 -5863.5625
[36,] -6296.3125 -4522.4375
[37,] -7045.6875 -6296.3125
[38,] -6467.3125 -7045.6875
[39,] -7575.8125 -6467.3125
[40,] -7514.6875 -7575.8125
[41,] -9361.9375 -7514.6875
[42,] -13668.8125 -9361.9375
[43,] -11217.3125 -13668.8125
[44,] -13481.3125 -11217.3125
[45,] -15594.0625 -13481.3125
[46,] -15611.5625 -15594.0625
[47,] -16458.4375 -15611.5625
[48,] -3886.1875 -16458.4375
[49,] -4538.5625 -3886.1875
[50,] -5664.1875 -4538.5625
[51,] -4542.6875 -5664.1875
[52,] -4958.5625 -4542.6875
[53,] -5550.8125 -4958.5625
[54,] -11003.6875 -5550.8125
[55,] -6340.1875 -11003.6875
[56,] -10047.1875 -6340.1875
[57,] -9435.9375 -10047.1875
[58,] -9628.4375 -9435.9375
[59,] -11148.3125 -9628.4375
[60,] -10702.1875 -11148.3125
[61,] -9817.5625 -10702.1875
[62,] -8089.1875 -9817.5625
[63,] -6839.6875 -8089.1875
[64,] -5167.5625 -6839.6875
[65,] -5664.8125 -5167.5625
[66,] -8067.6875 -5664.8125
[67,] -6418.1875 -8067.6875
[68,] -3899.1875 -6418.1875
[69,] -1526.9375 -3899.1875
[70,] -6261.4375 -1526.9375
[71,] -6858.3125 -6261.4375
[72,] -7096.1875 -6858.3125
[73,] -5010.5625 -7096.1875
[74,] -2412.1875 -5010.5625
[75,] -1684.6875 -2412.1875
[76,] -87.5625 -1684.6875
[77,] 1134.1875 -87.5625
[78,] 6952.3125 1134.1875
[79,] 1977.8125 6952.3125
[80,] 2784.8125 1977.8125
[81,] 5571.0625 2784.8125
[82,] 7273.5625 5571.0625
[83,] 7575.6875 7273.5625
[84,] 6950.8125 7575.6875
[85,] 7397.4375 6950.8125
[86,] 8284.8125 7397.4375
[87,] 9011.3125 8284.8125
[88,] 10105.4375 9011.3125
[89,] 14031.1875 10105.4375
[90,] 20414.3125 14031.1875
[91,] 13527.8125 20414.3125
[92,] 14141.8125 13527.8125
[93,] 14419.0625 14141.8125
[94,] 14440.5625 14419.0625
[95,] 16354.6875 14440.5625
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 12661.3125 14899.6875
2 9572.6875 12661.3125
3 9157.1875 9572.6875
4 7387.3125 9157.1875
5 6936.0625 7387.3125
6 10302.1875 6936.0625
7 9065.6875 10302.1875
8 8959.6875 9065.6875
9 9940.9375 8959.6875
10 11885.4375 9940.9375
11 11096.5625 11885.4375
12 6821.6875 11096.5625
13 6361.3125 6821.6875
14 5751.6875 6361.3125
15 5051.1875 5751.6875
16 3518.3125 5051.1875
17 4744.0625 3518.3125
18 2773.1875 4744.0625
19 -181.3125 2773.1875
20 4693.6875 -181.3125
21 3109.9375 4693.6875
22 3765.4375 3109.9375
23 3960.5625 3765.4375
24 -691.3125 3960.5625
25 -7.6875 -691.3125
26 -976.3125 -7.6875
27 -2576.8125 -976.3125
28 -3282.6875 -2576.8125
29 -6267.9375 -3282.6875
30 -7701.8125 -6267.9375
31 -414.3125 -7701.8125
32 -3152.3125 -414.3125
33 -6484.0625 -3152.3125
34 -5863.5625 -6484.0625
35 -4522.4375 -5863.5625
36 -6296.3125 -4522.4375
37 -7045.6875 -6296.3125
38 -6467.3125 -7045.6875
39 -7575.8125 -6467.3125
40 -7514.6875 -7575.8125
41 -9361.9375 -7514.6875
42 -13668.8125 -9361.9375
43 -11217.3125 -13668.8125
44 -13481.3125 -11217.3125
45 -15594.0625 -13481.3125
46 -15611.5625 -15594.0625
47 -16458.4375 -15611.5625
48 -3886.1875 -16458.4375
49 -4538.5625 -3886.1875
50 -5664.1875 -4538.5625
51 -4542.6875 -5664.1875
52 -4958.5625 -4542.6875
53 -5550.8125 -4958.5625
54 -11003.6875 -5550.8125
55 -6340.1875 -11003.6875
56 -10047.1875 -6340.1875
57 -9435.9375 -10047.1875
58 -9628.4375 -9435.9375
59 -11148.3125 -9628.4375
60 -10702.1875 -11148.3125
61 -9817.5625 -10702.1875
62 -8089.1875 -9817.5625
63 -6839.6875 -8089.1875
64 -5167.5625 -6839.6875
65 -5664.8125 -5167.5625
66 -8067.6875 -5664.8125
67 -6418.1875 -8067.6875
68 -3899.1875 -6418.1875
69 -1526.9375 -3899.1875
70 -6261.4375 -1526.9375
71 -6858.3125 -6261.4375
72 -7096.1875 -6858.3125
73 -5010.5625 -7096.1875
74 -2412.1875 -5010.5625
75 -1684.6875 -2412.1875
76 -87.5625 -1684.6875
77 1134.1875 -87.5625
78 6952.3125 1134.1875
79 1977.8125 6952.3125
80 2784.8125 1977.8125
81 5571.0625 2784.8125
82 7273.5625 5571.0625
83 7575.6875 7273.5625
84 6950.8125 7575.6875
85 7397.4375 6950.8125
86 8284.8125 7397.4375
87 9011.3125 8284.8125
88 10105.4375 9011.3125
89 14031.1875 10105.4375
90 20414.3125 14031.1875
91 13527.8125 20414.3125
92 14141.8125 13527.8125
93 14419.0625 14141.8125
94 14440.5625 14419.0625
95 16354.6875 14440.5625
> 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/7sdk01258545955.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/8eas51258545955.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/9gylc1258545955.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')
hat values (leverages) are all = 0.1354167
and there are no factor predictors; no plot no. 5
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/rcomp/tmp/10yexh1258545955.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/11w0k41258545955.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/12uqqq1258545955.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/13rp6n1258545955.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/142x3k1258545955.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/15807h1258545955.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/1605t81258545955.tab")
+ }
>
> system("convert tmp/1m7291258545954.ps tmp/1m7291258545954.png")
> system("convert tmp/2s6a41258545954.ps tmp/2s6a41258545954.png")
> system("convert tmp/3gl0n1258545954.ps tmp/3gl0n1258545954.png")
> system("convert tmp/45rm91258545954.ps tmp/45rm91258545954.png")
> system("convert tmp/5rxzx1258545955.ps tmp/5rxzx1258545955.png")
> system("convert tmp/60vup1258545955.ps tmp/60vup1258545955.png")
> system("convert tmp/7sdk01258545955.ps tmp/7sdk01258545955.png")
> system("convert tmp/8eas51258545955.ps tmp/8eas51258545955.png")
> system("convert tmp/9gylc1258545955.ps tmp/9gylc1258545955.png")
> system("convert tmp/10yexh1258545955.ps tmp/10yexh1258545955.png")
>
>
> proc.time()
user system elapsed
2.871 1.576 5.873