R version 2.8.0 (2008-10-20)
Copyright (C) 2008 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(98.5,0,96.7,0,113.1,0,100,0,104.7,0,108.5,0,90.5,0,88.6,0,105.4,0,119.9,0,107.2,0,84.1,0,101.4,0,105.1,0,118.7,0,113.8,0,113.8,0,118.9,0,98.5,0,91,0,120.7,0,127.9,0,112.4,0,93.1,0,107.5,0,107.3,0,114.8,0,120.8,0,112.2,0,123.3,0,100.6,0,86.7,0,123.6,0,125.3,0,111.1,0,98.4,0,102.3,0,105,0,128.2,0,124.7,0,116.1,0,131.2,0,97.7,0,88.8,0,132.8,0,113.9,0,112.6,1,104.3,1,107.5,1,106,1,117.3,1,123.1,1,114.3,1,132,1,92.3,1,93.7,1,121.3,1,113.6,1,116.3,1,98.3,1,111.9,1,109.3,1,133.2,1,118,1,131.6,1,134.1,1,96.7,1,99.8,1,128.3,1,134.9,1,130.7,1,107.3,1,121.6,1,120.6,1,140.5,1,124.8,1,129.9,1,159.4,1,111,1,110.1,1,132.7,1,135,1,118.6,1,94,1,117.9,1,114.7,1,113.6,1,130.6,1,117.1,1,123.2,1,106.1,1,87.9,1),dim=c(2,92),dimnames=list(c('y','x'),1:92))
> y <- array(NA,dim=c(2,92),dimnames=list(c('y','x'),1:92))
> 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 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 98.5 0 1 0 0 0 0 0 0 0 0 0 0 1
2 96.7 0 0 1 0 0 0 0 0 0 0 0 0 2
3 113.1 0 0 0 1 0 0 0 0 0 0 0 0 3
4 100.0 0 0 0 0 1 0 0 0 0 0 0 0 4
5 104.7 0 0 0 0 0 1 0 0 0 0 0 0 5
6 108.5 0 0 0 0 0 0 1 0 0 0 0 0 6
7 90.5 0 0 0 0 0 0 0 1 0 0 0 0 7
8 88.6 0 0 0 0 0 0 0 0 1 0 0 0 8
9 105.4 0 0 0 0 0 0 0 0 0 1 0 0 9
10 119.9 0 0 0 0 0 0 0 0 0 0 1 0 10
11 107.2 0 0 0 0 0 0 0 0 0 0 0 1 11
12 84.1 0 0 0 0 0 0 0 0 0 0 0 0 12
13 101.4 0 1 0 0 0 0 0 0 0 0 0 0 13
14 105.1 0 0 1 0 0 0 0 0 0 0 0 0 14
15 118.7 0 0 0 1 0 0 0 0 0 0 0 0 15
16 113.8 0 0 0 0 1 0 0 0 0 0 0 0 16
17 113.8 0 0 0 0 0 1 0 0 0 0 0 0 17
18 118.9 0 0 0 0 0 0 1 0 0 0 0 0 18
19 98.5 0 0 0 0 0 0 0 1 0 0 0 0 19
20 91.0 0 0 0 0 0 0 0 0 1 0 0 0 20
21 120.7 0 0 0 0 0 0 0 0 0 1 0 0 21
22 127.9 0 0 0 0 0 0 0 0 0 0 1 0 22
23 112.4 0 0 0 0 0 0 0 0 0 0 0 1 23
24 93.1 0 0 0 0 0 0 0 0 0 0 0 0 24
25 107.5 0 1 0 0 0 0 0 0 0 0 0 0 25
26 107.3 0 0 1 0 0 0 0 0 0 0 0 0 26
27 114.8 0 0 0 1 0 0 0 0 0 0 0 0 27
28 120.8 0 0 0 0 1 0 0 0 0 0 0 0 28
29 112.2 0 0 0 0 0 1 0 0 0 0 0 0 29
30 123.3 0 0 0 0 0 0 1 0 0 0 0 0 30
31 100.6 0 0 0 0 0 0 0 1 0 0 0 0 31
32 86.7 0 0 0 0 0 0 0 0 1 0 0 0 32
33 123.6 0 0 0 0 0 0 0 0 0 1 0 0 33
34 125.3 0 0 0 0 0 0 0 0 0 0 1 0 34
35 111.1 0 0 0 0 0 0 0 0 0 0 0 1 35
36 98.4 0 0 0 0 0 0 0 0 0 0 0 0 36
37 102.3 0 1 0 0 0 0 0 0 0 0 0 0 37
38 105.0 0 0 1 0 0 0 0 0 0 0 0 0 38
39 128.2 0 0 0 1 0 0 0 0 0 0 0 0 39
40 124.7 0 0 0 0 1 0 0 0 0 0 0 0 40
41 116.1 0 0 0 0 0 1 0 0 0 0 0 0 41
42 131.2 0 0 0 0 0 0 1 0 0 0 0 0 42
43 97.7 0 0 0 0 0 0 0 1 0 0 0 0 43
44 88.8 0 0 0 0 0 0 0 0 1 0 0 0 44
45 132.8 0 0 0 0 0 0 0 0 0 1 0 0 45
46 113.9 0 0 0 0 0 0 0 0 0 0 1 0 46
47 112.6 1 0 0 0 0 0 0 0 0 0 0 1 47
48 104.3 1 0 0 0 0 0 0 0 0 0 0 0 48
49 107.5 1 1 0 0 0 0 0 0 0 0 0 0 49
50 106.0 1 0 1 0 0 0 0 0 0 0 0 0 50
51 117.3 1 0 0 1 0 0 0 0 0 0 0 0 51
52 123.1 1 0 0 0 1 0 0 0 0 0 0 0 52
53 114.3 1 0 0 0 0 1 0 0 0 0 0 0 53
54 132.0 1 0 0 0 0 0 1 0 0 0 0 0 54
55 92.3 1 0 0 0 0 0 0 1 0 0 0 0 55
56 93.7 1 0 0 0 0 0 0 0 1 0 0 0 56
57 121.3 1 0 0 0 0 0 0 0 0 1 0 0 57
58 113.6 1 0 0 0 0 0 0 0 0 0 1 0 58
59 116.3 1 0 0 0 0 0 0 0 0 0 0 1 59
60 98.3 1 0 0 0 0 0 0 0 0 0 0 0 60
61 111.9 1 1 0 0 0 0 0 0 0 0 0 0 61
62 109.3 1 0 1 0 0 0 0 0 0 0 0 0 62
63 133.2 1 0 0 1 0 0 0 0 0 0 0 0 63
64 118.0 1 0 0 0 1 0 0 0 0 0 0 0 64
65 131.6 1 0 0 0 0 1 0 0 0 0 0 0 65
66 134.1 1 0 0 0 0 0 1 0 0 0 0 0 66
67 96.7 1 0 0 0 0 0 0 1 0 0 0 0 67
68 99.8 1 0 0 0 0 0 0 0 1 0 0 0 68
69 128.3 1 0 0 0 0 0 0 0 0 1 0 0 69
70 134.9 1 0 0 0 0 0 0 0 0 0 1 0 70
71 130.7 1 0 0 0 0 0 0 0 0 0 0 1 71
72 107.3 1 0 0 0 0 0 0 0 0 0 0 0 72
73 121.6 1 1 0 0 0 0 0 0 0 0 0 0 73
74 120.6 1 0 1 0 0 0 0 0 0 0 0 0 74
75 140.5 1 0 0 1 0 0 0 0 0 0 0 0 75
76 124.8 1 0 0 0 1 0 0 0 0 0 0 0 76
77 129.9 1 0 0 0 0 1 0 0 0 0 0 0 77
78 159.4 1 0 0 0 0 0 1 0 0 0 0 0 78
79 111.0 1 0 0 0 0 0 0 1 0 0 0 0 79
80 110.1 1 0 0 0 0 0 0 0 1 0 0 0 80
81 132.7 1 0 0 0 0 0 0 0 0 1 0 0 81
82 135.0 1 0 0 0 0 0 0 0 0 0 1 0 82
83 118.6 1 0 0 0 0 0 0 0 0 0 0 1 83
84 94.0 1 0 0 0 0 0 0 0 0 0 0 0 84
85 117.9 1 1 0 0 0 0 0 0 0 0 0 0 85
86 114.7 1 0 1 0 0 0 0 0 0 0 0 0 86
87 113.6 1 0 0 1 0 0 0 0 0 0 0 0 87
88 130.6 1 0 0 0 1 0 0 0 0 0 0 0 88
89 117.1 1 0 0 0 0 1 0 0 0 0 0 0 89
90 123.2 1 0 0 0 0 0 1 0 0 0 0 0 90
91 106.1 1 0 0 0 0 0 0 1 0 0 0 0 91
92 87.9 1 0 0 0 0 0 0 0 1 0 0 0 92
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) x M1 M2 M3 M4
86.4632 -1.8760 12.5863 11.8554 25.9496 22.7562
M5 M6 M7 M8 M9 M10
20.5004 31.6196 1.7262 -4.3671 26.9334 27.5044
M11 t
18.7291 0.2433
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-18.10721 -3.25725 -0.02689 3.40575 24.21284
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 86.46319 3.03081 28.528 < 2e-16 ***
x -1.87596 2.97079 -0.631 0.529580
M1 12.58626 3.63192 3.465 0.000864 ***
M2 11.85543 3.63089 3.265 0.001626 **
M3 25.94959 3.63071 7.147 4.13e-10 ***
M4 22.75625 3.63140 6.267 1.89e-08 ***
M5 20.50041 3.63295 5.643 2.59e-07 ***
M6 31.61957 3.63536 8.698 4.17e-13 ***
M7 1.72624 3.63863 0.474 0.636526
M8 -4.36710 3.64275 -1.199 0.234217
M9 26.93345 3.75943 7.164 3.83e-10 ***
M10 27.50440 3.76284 7.309 2.02e-10 ***
M11 18.72905 3.74860 4.996 3.50e-06 ***
t 0.24334 0.05592 4.352 4.05e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.012 on 78 degrees of freedom
Multiple R-squared: 0.7922, Adjusted R-squared: 0.7576
F-statistic: 22.88 on 13 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,] 1.093442e-01 0.2186883318 0.8906558
[2,] 4.590352e-02 0.0918070420 0.9540965
[3,] 1.576829e-02 0.0315365878 0.9842317
[4,] 1.282813e-02 0.0256562657 0.9871719
[5,] 1.553455e-02 0.0310691053 0.9844654
[6,] 6.637771e-03 0.0132755422 0.9933622
[7,] 3.096793e-03 0.0061935851 0.9969032
[8,] 1.146186e-03 0.0022923729 0.9988538
[9,] 7.495525e-04 0.0014991049 0.9992504
[10,] 5.257373e-04 0.0010514745 0.9994743
[11,] 3.358680e-03 0.0067173609 0.9966413
[12,] 2.300613e-03 0.0046012252 0.9976994
[13,] 2.044453e-03 0.0040889057 0.9979555
[14,] 1.004201e-03 0.0020084030 0.9989958
[15,] 5.349565e-04 0.0010699130 0.9994650
[16,] 1.836116e-03 0.0036722316 0.9981639
[17,] 9.889787e-04 0.0019779574 0.9990110
[18,] 8.662256e-04 0.0017324512 0.9991338
[19,] 6.734352e-04 0.0013468703 0.9993266
[20,] 3.583955e-04 0.0007167911 0.9996416
[21,] 6.445608e-04 0.0012891216 0.9993554
[22,] 5.326622e-04 0.0010653243 0.9994673
[23,] 3.817777e-04 0.0007635554 0.9996182
[24,] 2.709858e-04 0.0005419716 0.9997290
[25,] 1.425986e-04 0.0002851972 0.9998574
[26,] 1.026864e-04 0.0002053728 0.9998973
[27,] 1.003324e-04 0.0002006649 0.9998997
[28,] 1.072771e-04 0.0002145543 0.9998927
[29,] 1.795800e-04 0.0003591599 0.9998204
[30,] 1.678852e-03 0.0033577040 0.9983211
[31,] 9.784222e-04 0.0019568444 0.9990216
[32,] 9.248979e-04 0.0018497957 0.9990751
[33,] 5.789644e-04 0.0011579288 0.9994210
[34,] 3.804547e-04 0.0007609093 0.9996195
[35,] 3.494953e-04 0.0006989906 0.9996505
[36,] 1.951849e-04 0.0003903698 0.9998048
[37,] 1.312902e-04 0.0002625804 0.9998687
[38,] 9.313772e-05 0.0001862754 0.9999069
[39,] 1.399041e-04 0.0002798083 0.9998601
[40,] 7.354072e-05 0.0001470814 0.9999265
[41,] 5.103745e-05 0.0001020749 0.9999490
[42,] 3.698381e-04 0.0007396763 0.9996302
[43,] 2.765559e-04 0.0005531117 0.9997234
[44,] 1.517510e-04 0.0003035019 0.9998482
[45,] 1.240172e-04 0.0002480344 0.9998760
[46,] 1.174452e-04 0.0002348903 0.9998826
[47,] 8.966861e-05 0.0001793372 0.9999103
[48,] 1.466101e-04 0.0002932201 0.9998534
[49,] 1.844676e-04 0.0003689353 0.9998155
[50,] 2.682015e-04 0.0005364030 0.9997318
[51,] 2.200132e-03 0.0044002637 0.9977999
[52,] 3.227512e-03 0.0064550249 0.9967725
[53,] 4.956200e-03 0.0099124008 0.9950438
[54,] 6.649143e-03 0.0132982866 0.9933509
[55,] 4.689654e-03 0.0093793081 0.9953103
[56,] 2.202771e-03 0.0044055417 0.9977972
[57,] 1.870399e-03 0.0037407972 0.9981296
[58,] 1.311967e-03 0.0026239334 0.9986880
[59,] 1.309475e-03 0.0026189494 0.9986905
> postscript(file="/var/www/html/freestat/rcomp/tmp/1vepz1227731738.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/freestat/rcomp/tmp/2yduv1227731738.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/freestat/rcomp/tmp/38aub1227731738.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/freestat/rcomp/tmp/4c0g91227731738.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/freestat/rcomp/tmp/5pjr81227731738.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 = 92
Frequency = 1
1 2 3 4 5 6
-0.79278846 -2.10528846 -0.04278846 -10.19278846 -3.48028846 -11.04278846
7 8 9 10 11 12
0.60721154 4.55721154 -10.18667582 3.49903846 -0.66895604 -5.28324176
13 14 15 16 17 18
-0.81284341 3.37465659 2.63715659 0.68715659 2.69965659 -3.56284341
19 20 21 22 23 24
5.68715659 4.03715659 2.19326923 8.57898352 1.61098901 0.79670330
25 26 27 28 29 30
2.36710165 2.65460165 -4.18289835 4.76710165 -1.82039835 -2.08289835
31 32 33 34 35 36
4.86710165 -3.18289835 2.17321429 3.05892857 -2.60906593 3.17664835
37 38 39 40 41 42
-5.75295330 -2.56545330 6.29704670 5.74704670 -0.84045330 2.89704670
43 44 45 46 47 48
-0.95295330 -4.00295330 8.45315934 -11.26112637 -2.15315934 8.03255495
49 50 51 52 53 54
-1.59704670 -2.60954670 -5.64704670 3.10295330 -3.68454670 2.65295330
55 56 57 58 59 60
-7.39704670 -0.14704670 -4.09093407 -12.60521978 -1.37321429 -0.88750000
61 62 63 64 65 66
-0.11710165 -2.22960165 7.33289835 -4.91710165 10.69539835 1.83289835
67 68 69 70 71 72
-5.91710165 3.03289835 -0.01098901 5.77472527 10.10673077 5.19244505
73 74 75 76 77 78
6.66284341 6.15034341 11.71284341 -1.03715659 6.07534341 24.21284341
79 80 81 82 83 84
5.46284341 10.41284341 1.46895604 2.95467033 -4.91332418 -11.02760989
85 86 87 88 89 90
0.04278846 -2.66971154 -18.10721154 1.84278846 -9.64471154 -14.90721154
91 92
-2.35721154 -14.70721154
> postscript(file="/var/www/html/freestat/rcomp/tmp/6fctk1227731738.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 = 92
Frequency = 1
lag(myerror, k = 1) myerror
0 -0.79278846 NA
1 -2.10528846 -0.79278846
2 -0.04278846 -2.10528846
3 -10.19278846 -0.04278846
4 -3.48028846 -10.19278846
5 -11.04278846 -3.48028846
6 0.60721154 -11.04278846
7 4.55721154 0.60721154
8 -10.18667582 4.55721154
9 3.49903846 -10.18667582
10 -0.66895604 3.49903846
11 -5.28324176 -0.66895604
12 -0.81284341 -5.28324176
13 3.37465659 -0.81284341
14 2.63715659 3.37465659
15 0.68715659 2.63715659
16 2.69965659 0.68715659
17 -3.56284341 2.69965659
18 5.68715659 -3.56284341
19 4.03715659 5.68715659
20 2.19326923 4.03715659
21 8.57898352 2.19326923
22 1.61098901 8.57898352
23 0.79670330 1.61098901
24 2.36710165 0.79670330
25 2.65460165 2.36710165
26 -4.18289835 2.65460165
27 4.76710165 -4.18289835
28 -1.82039835 4.76710165
29 -2.08289835 -1.82039835
30 4.86710165 -2.08289835
31 -3.18289835 4.86710165
32 2.17321429 -3.18289835
33 3.05892857 2.17321429
34 -2.60906593 3.05892857
35 3.17664835 -2.60906593
36 -5.75295330 3.17664835
37 -2.56545330 -5.75295330
38 6.29704670 -2.56545330
39 5.74704670 6.29704670
40 -0.84045330 5.74704670
41 2.89704670 -0.84045330
42 -0.95295330 2.89704670
43 -4.00295330 -0.95295330
44 8.45315934 -4.00295330
45 -11.26112637 8.45315934
46 -2.15315934 -11.26112637
47 8.03255495 -2.15315934
48 -1.59704670 8.03255495
49 -2.60954670 -1.59704670
50 -5.64704670 -2.60954670
51 3.10295330 -5.64704670
52 -3.68454670 3.10295330
53 2.65295330 -3.68454670
54 -7.39704670 2.65295330
55 -0.14704670 -7.39704670
56 -4.09093407 -0.14704670
57 -12.60521978 -4.09093407
58 -1.37321429 -12.60521978
59 -0.88750000 -1.37321429
60 -0.11710165 -0.88750000
61 -2.22960165 -0.11710165
62 7.33289835 -2.22960165
63 -4.91710165 7.33289835
64 10.69539835 -4.91710165
65 1.83289835 10.69539835
66 -5.91710165 1.83289835
67 3.03289835 -5.91710165
68 -0.01098901 3.03289835
69 5.77472527 -0.01098901
70 10.10673077 5.77472527
71 5.19244505 10.10673077
72 6.66284341 5.19244505
73 6.15034341 6.66284341
74 11.71284341 6.15034341
75 -1.03715659 11.71284341
76 6.07534341 -1.03715659
77 24.21284341 6.07534341
78 5.46284341 24.21284341
79 10.41284341 5.46284341
80 1.46895604 10.41284341
81 2.95467033 1.46895604
82 -4.91332418 2.95467033
83 -11.02760989 -4.91332418
84 0.04278846 -11.02760989
85 -2.66971154 0.04278846
86 -18.10721154 -2.66971154
87 1.84278846 -18.10721154
88 -9.64471154 1.84278846
89 -14.90721154 -9.64471154
90 -2.35721154 -14.90721154
91 -14.70721154 -2.35721154
92 NA -14.70721154
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -2.10528846 -0.79278846
[2,] -0.04278846 -2.10528846
[3,] -10.19278846 -0.04278846
[4,] -3.48028846 -10.19278846
[5,] -11.04278846 -3.48028846
[6,] 0.60721154 -11.04278846
[7,] 4.55721154 0.60721154
[8,] -10.18667582 4.55721154
[9,] 3.49903846 -10.18667582
[10,] -0.66895604 3.49903846
[11,] -5.28324176 -0.66895604
[12,] -0.81284341 -5.28324176
[13,] 3.37465659 -0.81284341
[14,] 2.63715659 3.37465659
[15,] 0.68715659 2.63715659
[16,] 2.69965659 0.68715659
[17,] -3.56284341 2.69965659
[18,] 5.68715659 -3.56284341
[19,] 4.03715659 5.68715659
[20,] 2.19326923 4.03715659
[21,] 8.57898352 2.19326923
[22,] 1.61098901 8.57898352
[23,] 0.79670330 1.61098901
[24,] 2.36710165 0.79670330
[25,] 2.65460165 2.36710165
[26,] -4.18289835 2.65460165
[27,] 4.76710165 -4.18289835
[28,] -1.82039835 4.76710165
[29,] -2.08289835 -1.82039835
[30,] 4.86710165 -2.08289835
[31,] -3.18289835 4.86710165
[32,] 2.17321429 -3.18289835
[33,] 3.05892857 2.17321429
[34,] -2.60906593 3.05892857
[35,] 3.17664835 -2.60906593
[36,] -5.75295330 3.17664835
[37,] -2.56545330 -5.75295330
[38,] 6.29704670 -2.56545330
[39,] 5.74704670 6.29704670
[40,] -0.84045330 5.74704670
[41,] 2.89704670 -0.84045330
[42,] -0.95295330 2.89704670
[43,] -4.00295330 -0.95295330
[44,] 8.45315934 -4.00295330
[45,] -11.26112637 8.45315934
[46,] -2.15315934 -11.26112637
[47,] 8.03255495 -2.15315934
[48,] -1.59704670 8.03255495
[49,] -2.60954670 -1.59704670
[50,] -5.64704670 -2.60954670
[51,] 3.10295330 -5.64704670
[52,] -3.68454670 3.10295330
[53,] 2.65295330 -3.68454670
[54,] -7.39704670 2.65295330
[55,] -0.14704670 -7.39704670
[56,] -4.09093407 -0.14704670
[57,] -12.60521978 -4.09093407
[58,] -1.37321429 -12.60521978
[59,] -0.88750000 -1.37321429
[60,] -0.11710165 -0.88750000
[61,] -2.22960165 -0.11710165
[62,] 7.33289835 -2.22960165
[63,] -4.91710165 7.33289835
[64,] 10.69539835 -4.91710165
[65,] 1.83289835 10.69539835
[66,] -5.91710165 1.83289835
[67,] 3.03289835 -5.91710165
[68,] -0.01098901 3.03289835
[69,] 5.77472527 -0.01098901
[70,] 10.10673077 5.77472527
[71,] 5.19244505 10.10673077
[72,] 6.66284341 5.19244505
[73,] 6.15034341 6.66284341
[74,] 11.71284341 6.15034341
[75,] -1.03715659 11.71284341
[76,] 6.07534341 -1.03715659
[77,] 24.21284341 6.07534341
[78,] 5.46284341 24.21284341
[79,] 10.41284341 5.46284341
[80,] 1.46895604 10.41284341
[81,] 2.95467033 1.46895604
[82,] -4.91332418 2.95467033
[83,] -11.02760989 -4.91332418
[84,] 0.04278846 -11.02760989
[85,] -2.66971154 0.04278846
[86,] -18.10721154 -2.66971154
[87,] 1.84278846 -18.10721154
[88,] -9.64471154 1.84278846
[89,] -14.90721154 -9.64471154
[90,] -2.35721154 -14.90721154
[91,] -14.70721154 -2.35721154
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -2.10528846 -0.79278846
2 -0.04278846 -2.10528846
3 -10.19278846 -0.04278846
4 -3.48028846 -10.19278846
5 -11.04278846 -3.48028846
6 0.60721154 -11.04278846
7 4.55721154 0.60721154
8 -10.18667582 4.55721154
9 3.49903846 -10.18667582
10 -0.66895604 3.49903846
11 -5.28324176 -0.66895604
12 -0.81284341 -5.28324176
13 3.37465659 -0.81284341
14 2.63715659 3.37465659
15 0.68715659 2.63715659
16 2.69965659 0.68715659
17 -3.56284341 2.69965659
18 5.68715659 -3.56284341
19 4.03715659 5.68715659
20 2.19326923 4.03715659
21 8.57898352 2.19326923
22 1.61098901 8.57898352
23 0.79670330 1.61098901
24 2.36710165 0.79670330
25 2.65460165 2.36710165
26 -4.18289835 2.65460165
27 4.76710165 -4.18289835
28 -1.82039835 4.76710165
29 -2.08289835 -1.82039835
30 4.86710165 -2.08289835
31 -3.18289835 4.86710165
32 2.17321429 -3.18289835
33 3.05892857 2.17321429
34 -2.60906593 3.05892857
35 3.17664835 -2.60906593
36 -5.75295330 3.17664835
37 -2.56545330 -5.75295330
38 6.29704670 -2.56545330
39 5.74704670 6.29704670
40 -0.84045330 5.74704670
41 2.89704670 -0.84045330
42 -0.95295330 2.89704670
43 -4.00295330 -0.95295330
44 8.45315934 -4.00295330
45 -11.26112637 8.45315934
46 -2.15315934 -11.26112637
47 8.03255495 -2.15315934
48 -1.59704670 8.03255495
49 -2.60954670 -1.59704670
50 -5.64704670 -2.60954670
51 3.10295330 -5.64704670
52 -3.68454670 3.10295330
53 2.65295330 -3.68454670
54 -7.39704670 2.65295330
55 -0.14704670 -7.39704670
56 -4.09093407 -0.14704670
57 -12.60521978 -4.09093407
58 -1.37321429 -12.60521978
59 -0.88750000 -1.37321429
60 -0.11710165 -0.88750000
61 -2.22960165 -0.11710165
62 7.33289835 -2.22960165
63 -4.91710165 7.33289835
64 10.69539835 -4.91710165
65 1.83289835 10.69539835
66 -5.91710165 1.83289835
67 3.03289835 -5.91710165
68 -0.01098901 3.03289835
69 5.77472527 -0.01098901
70 10.10673077 5.77472527
71 5.19244505 10.10673077
72 6.66284341 5.19244505
73 6.15034341 6.66284341
74 11.71284341 6.15034341
75 -1.03715659 11.71284341
76 6.07534341 -1.03715659
77 24.21284341 6.07534341
78 5.46284341 24.21284341
79 10.41284341 5.46284341
80 1.46895604 10.41284341
81 2.95467033 1.46895604
82 -4.91332418 2.95467033
83 -11.02760989 -4.91332418
84 0.04278846 -11.02760989
85 -2.66971154 0.04278846
86 -18.10721154 -2.66971154
87 1.84278846 -18.10721154
88 -9.64471154 1.84278846
89 -14.90721154 -9.64471154
90 -2.35721154 -14.90721154
91 -14.70721154 -2.35721154
> 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/freestat/rcomp/tmp/7c00f1227731738.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/freestat/rcomp/tmp/85fc81227731738.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/freestat/rcomp/tmp/999bw1227731738.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/freestat/rcomp/tmp/1060tv1227731738.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/freestat/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/freestat/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/freestat/rcomp/tmp/11jgxk1227731738.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/freestat/rcomp/tmp/12rry61227731738.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/freestat/rcomp/tmp/130u721227731738.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/freestat/rcomp/tmp/14ufnj1227731738.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/freestat/rcomp/tmp/159byq1227731738.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/freestat/rcomp/tmp/16v8dk1227731739.tab")
+ }
>
> system("convert tmp/1vepz1227731738.ps tmp/1vepz1227731738.png")
> system("convert tmp/2yduv1227731738.ps tmp/2yduv1227731738.png")
> system("convert tmp/38aub1227731738.ps tmp/38aub1227731738.png")
> system("convert tmp/4c0g91227731738.ps tmp/4c0g91227731738.png")
> system("convert tmp/5pjr81227731738.ps tmp/5pjr81227731738.png")
> system("convert tmp/6fctk1227731738.ps tmp/6fctk1227731738.png")
> system("convert tmp/7c00f1227731738.ps tmp/7c00f1227731738.png")
> system("convert tmp/85fc81227731738.ps tmp/85fc81227731738.png")
> system("convert tmp/999bw1227731738.ps tmp/999bw1227731738.png")
> system("convert tmp/1060tv1227731738.ps tmp/1060tv1227731738.png")
>
>
> proc.time()
user system elapsed
4.089 2.509 4.460