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(28,30,32,31,32,34,32,35,29,33,33,34,32,27,26,22,22,23,23,19,15,5,1,11,18,20,21,19,20,19,18,16,16,16,15,11,10,6,1,8,10,9,6,8,14,4,13,13,16,18,16,15,13,19,15,17,17,13,12,13,13,16,17,14,8,8,8,9,5,11,10,14,18,17,14,15,13,17,17,17,17,21,20,11,18,20,18,21,21,20,18,17,17,18,11,15,13,16,16,12,10,8,6,8,10),dim=c(1,105),dimnames=list(c('X'),1:105))
> y <- array(NA,dim=c(1,105),dimnames=list(c('X'),1:105))
> 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
X M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 28 1 0 0 0 0 0 0 0 0 0 0 1
2 30 0 1 0 0 0 0 0 0 0 0 0 2
3 32 0 0 1 0 0 0 0 0 0 0 0 3
4 31 0 0 0 1 0 0 0 0 0 0 0 4
5 32 0 0 0 0 1 0 0 0 0 0 0 5
6 34 0 0 0 0 0 1 0 0 0 0 0 6
7 32 0 0 0 0 0 0 1 0 0 0 0 7
8 35 0 0 0 0 0 0 0 1 0 0 0 8
9 29 0 0 0 0 0 0 0 0 1 0 0 9
10 33 0 0 0 0 0 0 0 0 0 1 0 10
11 33 0 0 0 0 0 0 0 0 0 0 1 11
12 34 0 0 0 0 0 0 0 0 0 0 0 12
13 32 1 0 0 0 0 0 0 0 0 0 0 13
14 27 0 1 0 0 0 0 0 0 0 0 0 14
15 26 0 0 1 0 0 0 0 0 0 0 0 15
16 22 0 0 0 1 0 0 0 0 0 0 0 16
17 22 0 0 0 0 1 0 0 0 0 0 0 17
18 23 0 0 0 0 0 1 0 0 0 0 0 18
19 23 0 0 0 0 0 0 1 0 0 0 0 19
20 19 0 0 0 0 0 0 0 1 0 0 0 20
21 15 0 0 0 0 0 0 0 0 1 0 0 21
22 5 0 0 0 0 0 0 0 0 0 1 0 22
23 1 0 0 0 0 0 0 0 0 0 0 1 23
24 11 0 0 0 0 0 0 0 0 0 0 0 24
25 18 1 0 0 0 0 0 0 0 0 0 0 25
26 20 0 1 0 0 0 0 0 0 0 0 0 26
27 21 0 0 1 0 0 0 0 0 0 0 0 27
28 19 0 0 0 1 0 0 0 0 0 0 0 28
29 20 0 0 0 0 1 0 0 0 0 0 0 29
30 19 0 0 0 0 0 1 0 0 0 0 0 30
31 18 0 0 0 0 0 0 1 0 0 0 0 31
32 16 0 0 0 0 0 0 0 1 0 0 0 32
33 16 0 0 0 0 0 0 0 0 1 0 0 33
34 16 0 0 0 0 0 0 0 0 0 1 0 34
35 15 0 0 0 0 0 0 0 0 0 0 1 35
36 11 0 0 0 0 0 0 0 0 0 0 0 36
37 10 1 0 0 0 0 0 0 0 0 0 0 37
38 6 0 1 0 0 0 0 0 0 0 0 0 38
39 1 0 0 1 0 0 0 0 0 0 0 0 39
40 8 0 0 0 1 0 0 0 0 0 0 0 40
41 10 0 0 0 0 1 0 0 0 0 0 0 41
42 9 0 0 0 0 0 1 0 0 0 0 0 42
43 6 0 0 0 0 0 0 1 0 0 0 0 43
44 8 0 0 0 0 0 0 0 1 0 0 0 44
45 14 0 0 0 0 0 0 0 0 1 0 0 45
46 4 0 0 0 0 0 0 0 0 0 1 0 46
47 13 0 0 0 0 0 0 0 0 0 0 1 47
48 13 0 0 0 0 0 0 0 0 0 0 0 48
49 16 1 0 0 0 0 0 0 0 0 0 0 49
50 18 0 1 0 0 0 0 0 0 0 0 0 50
51 16 0 0 1 0 0 0 0 0 0 0 0 51
52 15 0 0 0 1 0 0 0 0 0 0 0 52
53 13 0 0 0 0 1 0 0 0 0 0 0 53
54 19 0 0 0 0 0 1 0 0 0 0 0 54
55 15 0 0 0 0 0 0 1 0 0 0 0 55
56 17 0 0 0 0 0 0 0 1 0 0 0 56
57 17 0 0 0 0 0 0 0 0 1 0 0 57
58 13 0 0 0 0 0 0 0 0 0 1 0 58
59 12 0 0 0 0 0 0 0 0 0 0 1 59
60 13 0 0 0 0 0 0 0 0 0 0 0 60
61 13 1 0 0 0 0 0 0 0 0 0 0 61
62 16 0 1 0 0 0 0 0 0 0 0 0 62
63 17 0 0 1 0 0 0 0 0 0 0 0 63
64 14 0 0 0 1 0 0 0 0 0 0 0 64
65 8 0 0 0 0 1 0 0 0 0 0 0 65
66 8 0 0 0 0 0 1 0 0 0 0 0 66
67 8 0 0 0 0 0 0 1 0 0 0 0 67
68 9 0 0 0 0 0 0 0 1 0 0 0 68
69 5 0 0 0 0 0 0 0 0 1 0 0 69
70 11 0 0 0 0 0 0 0 0 0 1 0 70
71 10 0 0 0 0 0 0 0 0 0 0 1 71
72 14 0 0 0 0 0 0 0 0 0 0 0 72
73 18 1 0 0 0 0 0 0 0 0 0 0 73
74 17 0 1 0 0 0 0 0 0 0 0 0 74
75 14 0 0 1 0 0 0 0 0 0 0 0 75
76 15 0 0 0 1 0 0 0 0 0 0 0 76
77 13 0 0 0 0 1 0 0 0 0 0 0 77
78 17 0 0 0 0 0 1 0 0 0 0 0 78
79 17 0 0 0 0 0 0 1 0 0 0 0 79
80 17 0 0 0 0 0 0 0 1 0 0 0 80
81 17 0 0 0 0 0 0 0 0 1 0 0 81
82 21 0 0 0 0 0 0 0 0 0 1 0 82
83 20 0 0 0 0 0 0 0 0 0 0 1 83
84 11 0 0 0 0 0 0 0 0 0 0 0 84
85 18 1 0 0 0 0 0 0 0 0 0 0 85
86 20 0 1 0 0 0 0 0 0 0 0 0 86
87 18 0 0 1 0 0 0 0 0 0 0 0 87
88 21 0 0 0 1 0 0 0 0 0 0 0 88
89 21 0 0 0 0 1 0 0 0 0 0 0 89
90 20 0 0 0 0 0 1 0 0 0 0 0 90
91 18 0 0 0 0 0 0 1 0 0 0 0 91
92 17 0 0 0 0 0 0 0 1 0 0 0 92
93 17 0 0 0 0 0 0 0 0 1 0 0 93
94 18 0 0 0 0 0 0 0 0 0 1 0 94
95 11 0 0 0 0 0 0 0 0 0 0 1 95
96 15 0 0 0 0 0 0 0 0 0 0 0 96
97 13 1 0 0 0 0 0 0 0 0 0 0 97
98 16 0 1 0 0 0 0 0 0 0 0 0 98
99 16 0 0 1 0 0 0 0 0 0 0 0 99
100 12 0 0 0 1 0 0 0 0 0 0 0 100
101 10 0 0 0 0 1 0 0 0 0 0 0 101
102 8 0 0 0 0 0 1 0 0 0 0 0 102
103 6 0 0 0 0 0 0 1 0 0 0 0 103
104 8 0 0 0 0 0 0 0 1 0 0 0 104
105 10 0 0 0 0 0 0 0 0 1 0 0 105
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) M1 M2 M3 M4 M5
21.9122 2.5776 3.1454 2.2688 1.9477 1.1822
M6 M7 M8 M9 M10 M11
2.1944 0.7623 1.2190 0.6757 -0.3717 -0.9984
t
-0.1234
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-18.3694 -2.9902 -0.5488 4.0330 13.5683
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 21.91216 2.72991 8.027 3.17e-12 ***
M1 2.57758 3.36677 0.766 0.446
M2 3.14540 3.36610 0.934 0.353
M3 2.26877 3.36558 0.674 0.502
M4 1.94770 3.36521 0.579 0.564
M5 1.18218 3.36499 0.351 0.726
M6 2.19444 3.36491 0.652 0.516
M7 0.76226 3.36499 0.227 0.821
M8 1.21897 3.36521 0.362 0.718
M9 0.67568 3.36558 0.201 0.841
M10 -0.37175 3.46275 -0.107 0.915
M11 -0.99837 3.46254 -0.288 0.774
t -0.12337 0.02236 -5.517 3.15e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.925 on 92 degrees of freedom
Multiple R-squared: 0.2725, Adjusted R-squared: 0.1776
F-statistic: 2.871 on 12 and 92 DF, p-value: 0.002116
> 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.3651498 0.730299577 0.6348502115
[2,] 0.3328528 0.665705610 0.6671471950
[3,] 0.3074173 0.614834584 0.6925827080
[4,] 0.2387471 0.477494263 0.7612528685
[5,] 0.3462366 0.692473135 0.6537634324
[6,] 0.3255086 0.651017122 0.6744914392
[7,] 0.8137118 0.372576354 0.1862881769
[8,] 0.9792194 0.041561291 0.0207806455
[9,] 0.9805734 0.038853242 0.0194266209
[10,] 0.9776661 0.044667808 0.0223339042
[11,] 0.9798689 0.040262132 0.0201310660
[12,] 0.9827832 0.034433604 0.0172168020
[13,] 0.9826000 0.034799925 0.0173999624
[14,] 0.9849765 0.030047093 0.0150235465
[15,] 0.9825526 0.034894843 0.0174474217
[16,] 0.9814546 0.037090845 0.0185454225
[17,] 0.9757848 0.048430350 0.0242151749
[18,] 0.9743308 0.051338418 0.0256692088
[19,] 0.9774190 0.045161952 0.0225809760
[20,] 0.9800341 0.039931797 0.0199658984
[21,] 0.9702298 0.059540324 0.0297701620
[22,] 0.9587937 0.082412574 0.0412062869
[23,] 0.9638815 0.072236996 0.0361184978
[24,] 0.9885350 0.022929985 0.0114649924
[25,] 0.9861471 0.027705895 0.0138529476
[26,] 0.9797566 0.040486842 0.0202434211
[27,] 0.9734465 0.053107034 0.0265535172
[28,] 0.9699216 0.060156869 0.0300784344
[29,] 0.9633149 0.073370202 0.0366851012
[30,] 0.9673990 0.065201923 0.0326009613
[31,] 0.9813518 0.037296325 0.0186481627
[32,] 0.9861321 0.027735734 0.0138678668
[33,] 0.9856856 0.028628763 0.0143143815
[34,] 0.9884930 0.023014049 0.0115070244
[35,] 0.9925143 0.014971478 0.0074857389
[36,] 0.9934047 0.013190629 0.0065953143
[37,] 0.9930213 0.013957483 0.0069787417
[38,] 0.9906644 0.018671103 0.0093355516
[39,] 0.9932851 0.013429831 0.0067149155
[40,] 0.9924544 0.015091280 0.0075456398
[41,] 0.9932212 0.013557508 0.0067787540
[42,] 0.9945185 0.010962949 0.0054814745
[43,] 0.9940033 0.011993428 0.0059967138
[44,] 0.9920200 0.015959941 0.0079799704
[45,] 0.9889311 0.022137804 0.0110689019
[46,] 0.9847520 0.030496038 0.0152480190
[47,] 0.9805226 0.038954784 0.0194773921
[48,] 0.9769606 0.046078796 0.0230393980
[49,] 0.9689548 0.062090375 0.0310451876
[50,] 0.9644864 0.071027261 0.0355136307
[51,] 0.9629458 0.074108315 0.0370541573
[52,] 0.9588767 0.082246530 0.0411232652
[53,] 0.9546961 0.090607777 0.0453038883
[54,] 0.9797350 0.040530032 0.0202650162
[55,] 0.9916847 0.016630622 0.0083153110
[56,] 0.9949472 0.010105540 0.0050527701
[57,] 0.9929777 0.014044534 0.0070222671
[58,] 0.9903582 0.019283570 0.0096417849
[59,] 0.9896102 0.020779541 0.0103897707
[60,] 0.9923193 0.015361396 0.0076806978
[61,] 0.9938996 0.012200771 0.0061003857
[62,] 0.9973872 0.005225623 0.0026128113
[63,] 0.9966770 0.006646048 0.0033230241
[64,] 0.9948328 0.010334367 0.0051671835
[65,] 0.9927897 0.014420693 0.0072103466
[66,] 0.9932765 0.013446968 0.0067234840
[67,] 0.9911653 0.017669411 0.0088347054
[68,] 0.9851212 0.029757623 0.0148788113
[69,] 0.9977885 0.004423071 0.0022115353
[70,] 0.9953734 0.009253231 0.0046266157
[71,] 0.9934298 0.013140483 0.0065702414
[72,] 0.9991303 0.001739353 0.0008696763
[73,] 0.9960809 0.007838146 0.0039190728
[74,] 0.9821741 0.035651783 0.0178258914
> postscript(file="/var/www/html/freestat/rcomp/tmp/1ewrq1227445610.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/2pfzl1227445610.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/3kf4u1227445611.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/4j5xa1227445611.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/5vyao1227445611.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 = 105
Frequency = 1
1 2 3 4 5 6
3.6336336 5.1891892 8.1891892 7.6336336 9.5225225 10.6336336
7 8 9 10 11 12
10.1891892 12.8558559 7.5225225 12.6933183 13.4433183 13.5683183
13 14 15 16 17 18
9.1141141 3.6696697 3.6696697 0.1141141 1.0030030 1.1141141
19 20 21 22 23 24
2.6696697 -1.6636637 -4.9969970 -13.8262012 -17.0762012 -7.9512012
25 26 27 28 29 30
-3.4054054 -1.8498498 0.1501502 -1.4054054 0.4834835 -1.4054054
31 32 33 34 35 36
-0.8498498 -3.1831832 -2.5165165 -1.3457207 -1.5957207 -6.4707207
37 38 39 40 41 42
-9.9249249 -14.3693694 -18.3693694 -10.9249249 -8.0360360 -9.9249249
43 44 45 46 47 48
-11.3693694 -9.7027027 -3.0360360 -11.8652402 -2.1152402 -2.9902402
49 50 51 52 53 54
-2.4444444 -0.8888889 -1.8888889 -2.4444444 -3.5555556 1.5555556
55 56 57 58 59 60
-0.8888889 0.7777778 1.4444444 -1.3847598 -1.6347598 -1.5097598
61 62 63 64 65 66
-3.9639640 -1.4084084 0.5915916 -1.9639640 -7.0750751 -7.9639640
67 68 69 70 71 72
-6.4084084 -5.7417417 -9.0750751 -1.9042793 -2.1542793 0.9707207
73 74 75 76 77 78
2.5165165 1.0720721 -0.9279279 0.5165165 -0.5945946 2.5165165
79 80 81 82 83 84
4.0720721 3.7387387 4.4054054 9.5762012 9.3262012 -0.5487988
85 86 87 88 89 90
3.9969970 5.5525526 4.5525526 7.9969970 8.8858859 6.9969970
91 92 93 94 95 96
6.5525526 5.2192192 5.8858859 8.0566817 1.8066817 4.9316817
97 98 99 100 101 102
0.4774775 3.0330330 4.0330330 0.4774775 -0.6336336 -3.5225225
103 104 105
-3.9669670 -2.3003003 0.3663664
> postscript(file="/var/www/html/freestat/rcomp/tmp/6hkrf1227445611.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 = 105
Frequency = 1
lag(myerror, k = 1) myerror
0 3.6336336 NA
1 5.1891892 3.6336336
2 8.1891892 5.1891892
3 7.6336336 8.1891892
4 9.5225225 7.6336336
5 10.6336336 9.5225225
6 10.1891892 10.6336336
7 12.8558559 10.1891892
8 7.5225225 12.8558559
9 12.6933183 7.5225225
10 13.4433183 12.6933183
11 13.5683183 13.4433183
12 9.1141141 13.5683183
13 3.6696697 9.1141141
14 3.6696697 3.6696697
15 0.1141141 3.6696697
16 1.0030030 0.1141141
17 1.1141141 1.0030030
18 2.6696697 1.1141141
19 -1.6636637 2.6696697
20 -4.9969970 -1.6636637
21 -13.8262012 -4.9969970
22 -17.0762012 -13.8262012
23 -7.9512012 -17.0762012
24 -3.4054054 -7.9512012
25 -1.8498498 -3.4054054
26 0.1501502 -1.8498498
27 -1.4054054 0.1501502
28 0.4834835 -1.4054054
29 -1.4054054 0.4834835
30 -0.8498498 -1.4054054
31 -3.1831832 -0.8498498
32 -2.5165165 -3.1831832
33 -1.3457207 -2.5165165
34 -1.5957207 -1.3457207
35 -6.4707207 -1.5957207
36 -9.9249249 -6.4707207
37 -14.3693694 -9.9249249
38 -18.3693694 -14.3693694
39 -10.9249249 -18.3693694
40 -8.0360360 -10.9249249
41 -9.9249249 -8.0360360
42 -11.3693694 -9.9249249
43 -9.7027027 -11.3693694
44 -3.0360360 -9.7027027
45 -11.8652402 -3.0360360
46 -2.1152402 -11.8652402
47 -2.9902402 -2.1152402
48 -2.4444444 -2.9902402
49 -0.8888889 -2.4444444
50 -1.8888889 -0.8888889
51 -2.4444444 -1.8888889
52 -3.5555556 -2.4444444
53 1.5555556 -3.5555556
54 -0.8888889 1.5555556
55 0.7777778 -0.8888889
56 1.4444444 0.7777778
57 -1.3847598 1.4444444
58 -1.6347598 -1.3847598
59 -1.5097598 -1.6347598
60 -3.9639640 -1.5097598
61 -1.4084084 -3.9639640
62 0.5915916 -1.4084084
63 -1.9639640 0.5915916
64 -7.0750751 -1.9639640
65 -7.9639640 -7.0750751
66 -6.4084084 -7.9639640
67 -5.7417417 -6.4084084
68 -9.0750751 -5.7417417
69 -1.9042793 -9.0750751
70 -2.1542793 -1.9042793
71 0.9707207 -2.1542793
72 2.5165165 0.9707207
73 1.0720721 2.5165165
74 -0.9279279 1.0720721
75 0.5165165 -0.9279279
76 -0.5945946 0.5165165
77 2.5165165 -0.5945946
78 4.0720721 2.5165165
79 3.7387387 4.0720721
80 4.4054054 3.7387387
81 9.5762012 4.4054054
82 9.3262012 9.5762012
83 -0.5487988 9.3262012
84 3.9969970 -0.5487988
85 5.5525526 3.9969970
86 4.5525526 5.5525526
87 7.9969970 4.5525526
88 8.8858859 7.9969970
89 6.9969970 8.8858859
90 6.5525526 6.9969970
91 5.2192192 6.5525526
92 5.8858859 5.2192192
93 8.0566817 5.8858859
94 1.8066817 8.0566817
95 4.9316817 1.8066817
96 0.4774775 4.9316817
97 3.0330330 0.4774775
98 4.0330330 3.0330330
99 0.4774775 4.0330330
100 -0.6336336 0.4774775
101 -3.5225225 -0.6336336
102 -3.9669670 -3.5225225
103 -2.3003003 -3.9669670
104 0.3663664 -2.3003003
105 NA 0.3663664
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 5.1891892 3.6336336
[2,] 8.1891892 5.1891892
[3,] 7.6336336 8.1891892
[4,] 9.5225225 7.6336336
[5,] 10.6336336 9.5225225
[6,] 10.1891892 10.6336336
[7,] 12.8558559 10.1891892
[8,] 7.5225225 12.8558559
[9,] 12.6933183 7.5225225
[10,] 13.4433183 12.6933183
[11,] 13.5683183 13.4433183
[12,] 9.1141141 13.5683183
[13,] 3.6696697 9.1141141
[14,] 3.6696697 3.6696697
[15,] 0.1141141 3.6696697
[16,] 1.0030030 0.1141141
[17,] 1.1141141 1.0030030
[18,] 2.6696697 1.1141141
[19,] -1.6636637 2.6696697
[20,] -4.9969970 -1.6636637
[21,] -13.8262012 -4.9969970
[22,] -17.0762012 -13.8262012
[23,] -7.9512012 -17.0762012
[24,] -3.4054054 -7.9512012
[25,] -1.8498498 -3.4054054
[26,] 0.1501502 -1.8498498
[27,] -1.4054054 0.1501502
[28,] 0.4834835 -1.4054054
[29,] -1.4054054 0.4834835
[30,] -0.8498498 -1.4054054
[31,] -3.1831832 -0.8498498
[32,] -2.5165165 -3.1831832
[33,] -1.3457207 -2.5165165
[34,] -1.5957207 -1.3457207
[35,] -6.4707207 -1.5957207
[36,] -9.9249249 -6.4707207
[37,] -14.3693694 -9.9249249
[38,] -18.3693694 -14.3693694
[39,] -10.9249249 -18.3693694
[40,] -8.0360360 -10.9249249
[41,] -9.9249249 -8.0360360
[42,] -11.3693694 -9.9249249
[43,] -9.7027027 -11.3693694
[44,] -3.0360360 -9.7027027
[45,] -11.8652402 -3.0360360
[46,] -2.1152402 -11.8652402
[47,] -2.9902402 -2.1152402
[48,] -2.4444444 -2.9902402
[49,] -0.8888889 -2.4444444
[50,] -1.8888889 -0.8888889
[51,] -2.4444444 -1.8888889
[52,] -3.5555556 -2.4444444
[53,] 1.5555556 -3.5555556
[54,] -0.8888889 1.5555556
[55,] 0.7777778 -0.8888889
[56,] 1.4444444 0.7777778
[57,] -1.3847598 1.4444444
[58,] -1.6347598 -1.3847598
[59,] -1.5097598 -1.6347598
[60,] -3.9639640 -1.5097598
[61,] -1.4084084 -3.9639640
[62,] 0.5915916 -1.4084084
[63,] -1.9639640 0.5915916
[64,] -7.0750751 -1.9639640
[65,] -7.9639640 -7.0750751
[66,] -6.4084084 -7.9639640
[67,] -5.7417417 -6.4084084
[68,] -9.0750751 -5.7417417
[69,] -1.9042793 -9.0750751
[70,] -2.1542793 -1.9042793
[71,] 0.9707207 -2.1542793
[72,] 2.5165165 0.9707207
[73,] 1.0720721 2.5165165
[74,] -0.9279279 1.0720721
[75,] 0.5165165 -0.9279279
[76,] -0.5945946 0.5165165
[77,] 2.5165165 -0.5945946
[78,] 4.0720721 2.5165165
[79,] 3.7387387 4.0720721
[80,] 4.4054054 3.7387387
[81,] 9.5762012 4.4054054
[82,] 9.3262012 9.5762012
[83,] -0.5487988 9.3262012
[84,] 3.9969970 -0.5487988
[85,] 5.5525526 3.9969970
[86,] 4.5525526 5.5525526
[87,] 7.9969970 4.5525526
[88,] 8.8858859 7.9969970
[89,] 6.9969970 8.8858859
[90,] 6.5525526 6.9969970
[91,] 5.2192192 6.5525526
[92,] 5.8858859 5.2192192
[93,] 8.0566817 5.8858859
[94,] 1.8066817 8.0566817
[95,] 4.9316817 1.8066817
[96,] 0.4774775 4.9316817
[97,] 3.0330330 0.4774775
[98,] 4.0330330 3.0330330
[99,] 0.4774775 4.0330330
[100,] -0.6336336 0.4774775
[101,] -3.5225225 -0.6336336
[102,] -3.9669670 -3.5225225
[103,] -2.3003003 -3.9669670
[104,] 0.3663664 -2.3003003
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 5.1891892 3.6336336
2 8.1891892 5.1891892
3 7.6336336 8.1891892
4 9.5225225 7.6336336
5 10.6336336 9.5225225
6 10.1891892 10.6336336
7 12.8558559 10.1891892
8 7.5225225 12.8558559
9 12.6933183 7.5225225
10 13.4433183 12.6933183
11 13.5683183 13.4433183
12 9.1141141 13.5683183
13 3.6696697 9.1141141
14 3.6696697 3.6696697
15 0.1141141 3.6696697
16 1.0030030 0.1141141
17 1.1141141 1.0030030
18 2.6696697 1.1141141
19 -1.6636637 2.6696697
20 -4.9969970 -1.6636637
21 -13.8262012 -4.9969970
22 -17.0762012 -13.8262012
23 -7.9512012 -17.0762012
24 -3.4054054 -7.9512012
25 -1.8498498 -3.4054054
26 0.1501502 -1.8498498
27 -1.4054054 0.1501502
28 0.4834835 -1.4054054
29 -1.4054054 0.4834835
30 -0.8498498 -1.4054054
31 -3.1831832 -0.8498498
32 -2.5165165 -3.1831832
33 -1.3457207 -2.5165165
34 -1.5957207 -1.3457207
35 -6.4707207 -1.5957207
36 -9.9249249 -6.4707207
37 -14.3693694 -9.9249249
38 -18.3693694 -14.3693694
39 -10.9249249 -18.3693694
40 -8.0360360 -10.9249249
41 -9.9249249 -8.0360360
42 -11.3693694 -9.9249249
43 -9.7027027 -11.3693694
44 -3.0360360 -9.7027027
45 -11.8652402 -3.0360360
46 -2.1152402 -11.8652402
47 -2.9902402 -2.1152402
48 -2.4444444 -2.9902402
49 -0.8888889 -2.4444444
50 -1.8888889 -0.8888889
51 -2.4444444 -1.8888889
52 -3.5555556 -2.4444444
53 1.5555556 -3.5555556
54 -0.8888889 1.5555556
55 0.7777778 -0.8888889
56 1.4444444 0.7777778
57 -1.3847598 1.4444444
58 -1.6347598 -1.3847598
59 -1.5097598 -1.6347598
60 -3.9639640 -1.5097598
61 -1.4084084 -3.9639640
62 0.5915916 -1.4084084
63 -1.9639640 0.5915916
64 -7.0750751 -1.9639640
65 -7.9639640 -7.0750751
66 -6.4084084 -7.9639640
67 -5.7417417 -6.4084084
68 -9.0750751 -5.7417417
69 -1.9042793 -9.0750751
70 -2.1542793 -1.9042793
71 0.9707207 -2.1542793
72 2.5165165 0.9707207
73 1.0720721 2.5165165
74 -0.9279279 1.0720721
75 0.5165165 -0.9279279
76 -0.5945946 0.5165165
77 2.5165165 -0.5945946
78 4.0720721 2.5165165
79 3.7387387 4.0720721
80 4.4054054 3.7387387
81 9.5762012 4.4054054
82 9.3262012 9.5762012
83 -0.5487988 9.3262012
84 3.9969970 -0.5487988
85 5.5525526 3.9969970
86 4.5525526 5.5525526
87 7.9969970 4.5525526
88 8.8858859 7.9969970
89 6.9969970 8.8858859
90 6.5525526 6.9969970
91 5.2192192 6.5525526
92 5.8858859 5.2192192
93 8.0566817 5.8858859
94 1.8066817 8.0566817
95 4.9316817 1.8066817
96 0.4774775 4.9316817
97 3.0330330 0.4774775
98 4.0330330 3.0330330
99 0.4774775 4.0330330
100 -0.6336336 0.4774775
101 -3.5225225 -0.6336336
102 -3.9669670 -3.5225225
103 -2.3003003 -3.9669670
104 0.3663664 -2.3003003
> 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/789kr1227445611.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/8mdp21227445611.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/97iyy1227445611.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/10ekwc1227445611.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/11t2cx1227445611.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/12731n1227445611.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/133q5c1227445611.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/14pmtd1227445611.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/15vd2s1227445611.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/163bpg1227445611.tab")
+ }
>
> system("convert tmp/1ewrq1227445610.ps tmp/1ewrq1227445610.png")
> system("convert tmp/2pfzl1227445610.ps tmp/2pfzl1227445610.png")
> system("convert tmp/3kf4u1227445611.ps tmp/3kf4u1227445611.png")
> system("convert tmp/4j5xa1227445611.ps tmp/4j5xa1227445611.png")
> system("convert tmp/5vyao1227445611.ps tmp/5vyao1227445611.png")
> system("convert tmp/6hkrf1227445611.ps tmp/6hkrf1227445611.png")
> system("convert tmp/789kr1227445611.ps tmp/789kr1227445611.png")
> system("convert tmp/8mdp21227445611.ps tmp/8mdp21227445611.png")
> system("convert tmp/97iyy1227445611.ps tmp/97iyy1227445611.png")
> system("convert tmp/10ekwc1227445611.ps tmp/10ekwc1227445611.png")
>
>
> proc.time()
user system elapsed
4.389 2.595 4.734