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.
Natural language support but running in an English locale
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(1
+ ,162556
+ ,807
+ ,213118
+ ,1
+ ,29790
+ ,444
+ ,81767
+ ,1
+ ,87550
+ ,412
+ ,153198
+ ,0
+ ,84738
+ ,428
+ ,-26007
+ ,1
+ ,54660
+ ,315
+ ,126942
+ ,1
+ ,42634
+ ,168
+ ,157214
+ ,0
+ ,40949
+ ,263
+ ,129352
+ ,1
+ ,45187
+ ,267
+ ,234817
+ ,1
+ ,37704
+ ,228
+ ,60448
+ ,1
+ ,16275
+ ,129
+ ,47818
+ ,0
+ ,25830
+ ,104
+ ,245546
+ ,0
+ ,12679
+ ,122
+ ,48020
+ ,1
+ ,18014
+ ,393
+ ,-1710
+ ,0
+ ,43556
+ ,190
+ ,32648
+ ,1
+ ,24811
+ ,280
+ ,95350
+ ,0
+ ,6575
+ ,63
+ ,151352
+ ,0
+ ,7123
+ ,102
+ ,288170
+ ,1
+ ,21950
+ ,265
+ ,114337
+ ,1
+ ,37597
+ ,234
+ ,37884
+ ,0
+ ,17821
+ ,277
+ ,122844
+ ,1
+ ,12988
+ ,73
+ ,82340
+ ,1
+ ,22330
+ ,67
+ ,79801
+ ,0
+ ,13326
+ ,103
+ ,165548
+ ,0
+ ,16189
+ ,290
+ ,116384
+ ,0
+ ,7146
+ ,83
+ ,134028
+ ,0
+ ,15824
+ ,56
+ ,63838
+ ,1
+ ,27664
+ ,236
+ ,74996
+ ,0
+ ,11920
+ ,73
+ ,31080
+ ,0
+ ,8568
+ ,34
+ ,32168
+ ,0
+ ,14416
+ ,139
+ ,49857
+ ,1
+ ,3369
+ ,26
+ ,87161
+ ,1
+ ,11819
+ ,70
+ ,106113
+ ,1
+ ,6984
+ ,40
+ ,80570
+ ,1
+ ,4519
+ ,42
+ ,102129
+ ,0
+ ,2220
+ ,12
+ ,301670
+ ,0
+ ,18562
+ ,211
+ ,102313
+ ,0
+ ,10327
+ ,74
+ ,88577
+ ,1
+ ,5336
+ ,80
+ ,112477
+ ,1
+ ,2365
+ ,83
+ ,191778
+ ,0
+ ,4069
+ ,131
+ ,79804
+ ,0
+ ,8636
+ ,203
+ ,128294
+ ,0
+ ,13718
+ ,56
+ ,96448
+ ,0
+ ,4525
+ ,89
+ ,93811
+ ,0
+ ,6869
+ ,88
+ ,117520
+ ,0
+ ,4628
+ ,39
+ ,69159
+ ,1
+ ,3689
+ ,25
+ ,101792
+ ,1
+ ,4891
+ ,49
+ ,210568
+ ,1
+ ,7489
+ ,149
+ ,136996
+ ,0
+ ,4901
+ ,58
+ ,121920
+ ,0
+ ,2284
+ ,41
+ ,76403
+ ,1
+ ,3160
+ ,90
+ ,108094
+ ,1
+ ,4150
+ ,136
+ ,134759
+ ,1
+ ,7285
+ ,97
+ ,188873
+ ,1
+ ,1134
+ ,63
+ ,146216
+ ,1
+ ,4658
+ ,114
+ ,156608
+ ,0
+ ,2384
+ ,77
+ ,61348
+ ,0
+ ,3748
+ ,6
+ ,50350
+ ,0
+ ,5371
+ ,47
+ ,87720
+ ,0
+ ,1285
+ ,51
+ ,99489
+ ,1
+ ,9327
+ ,85
+ ,87419
+ ,1
+ ,5565
+ ,43
+ ,94355
+ ,0
+ ,1528
+ ,32
+ ,60326
+ ,1
+ ,3122
+ ,25
+ ,94670
+ ,1
+ ,7561
+ ,77
+ ,82425
+ ,0
+ ,2675
+ ,54
+ ,59017
+ ,0
+ ,13253
+ ,251
+ ,90829
+ ,0
+ ,880
+ ,15
+ ,80791
+ ,1
+ ,2053
+ ,44
+ ,100423
+ ,0
+ ,1424
+ ,73
+ ,131116
+ ,1
+ ,4036
+ ,85
+ ,100269
+ ,1
+ ,3045
+ ,49
+ ,27330
+ ,0
+ ,5119
+ ,38
+ ,39039
+ ,0
+ ,1431
+ ,35
+ ,106885
+ ,0
+ ,554
+ ,9
+ ,79285
+ ,0
+ ,1975
+ ,34
+ ,118881
+ ,1
+ ,1765
+ ,20
+ ,77623
+ ,0
+ ,1012
+ ,29
+ ,114768
+ ,0
+ ,810
+ ,11
+ ,74015
+ ,0
+ ,1280
+ ,52
+ ,69465
+ ,1
+ ,666
+ ,13
+ ,117869
+ ,0
+ ,1380
+ ,29
+ ,60982
+ ,1
+ ,4677
+ ,66
+ ,90131
+ ,0
+ ,876
+ ,33
+ ,138971
+ ,0
+ ,814
+ ,15
+ ,39625
+ ,0
+ ,514
+ ,15
+ ,102725
+ ,1
+ ,5692
+ ,68
+ ,64239
+ ,0
+ ,3642
+ ,100
+ ,90262
+ ,0
+ ,540
+ ,13
+ ,103960
+ ,0
+ ,2099
+ ,45
+ ,106611
+ ,0
+ ,567
+ ,14
+ ,103345
+ ,0
+ ,2001
+ ,36
+ ,95551
+ ,1
+ ,2949
+ ,40
+ ,82903
+ ,0
+ ,2253
+ ,68
+ ,63593
+ ,1
+ ,6533
+ ,29
+ ,126910
+ ,0
+ ,1889
+ ,43
+ ,37527
+ ,1
+ ,3055
+ ,30
+ ,60247
+ ,0
+ ,272
+ ,9
+ ,112995
+ ,1
+ ,1414
+ ,22
+ ,70184
+ ,0
+ ,2564
+ ,19
+ ,130140
+ ,1
+ ,1383
+ ,9
+ ,73221)
+ ,dim=c(4
+ ,100)
+ ,dimnames=list(c('Group'
+ ,'Costs'
+ ,'Orders'
+ ,'Dividends')
+ ,1:100))
> y <- array(NA,dim=c(4,100),dimnames=list(c('Group','Costs','Orders','Dividends'),1:100))
> 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 = '3'
> #'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
Orders Group Costs Dividends
1 807 1 162556 213118
2 444 1 29790 81767
3 412 1 87550 153198
4 428 0 84738 -26007
5 315 1 54660 126942
6 168 1 42634 157214
7 263 0 40949 129352
8 267 1 45187 234817
9 228 1 37704 60448
10 129 1 16275 47818
11 104 0 25830 245546
12 122 0 12679 48020
13 393 1 18014 -1710
14 190 0 43556 32648
15 280 1 24811 95350
16 63 0 6575 151352
17 102 0 7123 288170
18 265 1 21950 114337
19 234 1 37597 37884
20 277 0 17821 122844
21 73 1 12988 82340
22 67 1 22330 79801
23 103 0 13326 165548
24 290 0 16189 116384
25 83 0 7146 134028
26 56 0 15824 63838
27 236 1 27664 74996
28 73 0 11920 31080
29 34 0 8568 32168
30 139 0 14416 49857
31 26 1 3369 87161
32 70 1 11819 106113
33 40 1 6984 80570
34 42 1 4519 102129
35 12 0 2220 301670
36 211 0 18562 102313
37 74 0 10327 88577
38 80 1 5336 112477
39 83 1 2365 191778
40 131 0 4069 79804
41 203 0 8636 128294
42 56 0 13718 96448
43 89 0 4525 93811
44 88 0 6869 117520
45 39 0 4628 69159
46 25 1 3689 101792
47 49 1 4891 210568
48 149 1 7489 136996
49 58 0 4901 121920
50 41 0 2284 76403
51 90 1 3160 108094
52 136 1 4150 134759
53 97 1 7285 188873
54 63 1 1134 146216
55 114 1 4658 156608
56 77 0 2384 61348
57 6 0 3748 50350
58 47 0 5371 87720
59 51 0 1285 99489
60 85 1 9327 87419
61 43 1 5565 94355
62 32 0 1528 60326
63 25 1 3122 94670
64 77 1 7561 82425
65 54 0 2675 59017
66 251 0 13253 90829
67 15 0 880 80791
68 44 1 2053 100423
69 73 0 1424 131116
70 85 1 4036 100269
71 49 1 3045 27330
72 38 0 5119 39039
73 35 0 1431 106885
74 9 0 554 79285
75 34 0 1975 118881
76 20 1 1765 77623
77 29 0 1012 114768
78 11 0 810 74015
79 52 0 1280 69465
80 13 1 666 117869
81 29 0 1380 60982
82 66 1 4677 90131
83 33 0 876 138971
84 15 0 814 39625
85 15 0 514 102725
86 68 1 5692 64239
87 100 0 3642 90262
88 13 0 540 103960
89 45 0 2099 106611
90 14 0 567 103345
91 36 0 2001 95551
92 40 1 2949 82903
93 68 0 2253 63593
94 29 1 6533 126910
95 43 0 1889 37527
96 30 1 3055 60247
97 9 0 272 112995
98 22 1 1414 70184
99 19 0 2564 130140
100 9 1 1383 73221
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Group Costs Dividends
4.579e+01 1.010e+01 4.878e-03 -4.408e-05
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-94.29 -33.25 -12.94 13.93 249.16
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.579e+01 1.385e+01 3.306 0.00133 **
Group 1.010e+01 1.242e+01 0.813 0.41821
Costs 4.878e-03 2.856e-04 17.081 < 2e-16 ***
Dividends -4.408e-05 1.149e-04 -0.384 0.70207
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 60.48 on 96 degrees of freedom
Multiple R-squared: 0.7651, Adjusted R-squared: 0.7578
F-statistic: 104.2 on 3 and 96 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.9993803 1.239397e-03 6.196986e-04
[2,] 0.9983580 3.283984e-03 1.641992e-03
[3,] 0.9981051 3.789855e-03 1.894927e-03
[4,] 0.9972015 5.596974e-03 2.798487e-03
[5,] 0.9962515 7.497003e-03 3.748502e-03
[6,] 0.9925740 1.485200e-02 7.426002e-03
[7,] 0.9999549 9.013038e-05 4.506519e-05
[8,] 0.9999877 2.452518e-05 1.226259e-05
[9,] 0.9999841 3.171414e-05 1.585707e-05
[10,] 0.9999668 6.641407e-05 3.320703e-05
[11,] 0.9999720 5.590913e-05 2.795456e-05
[12,] 0.9999693 6.133100e-05 3.066550e-05
[13,] 0.9999753 4.932030e-05 2.466015e-05
[14,] 0.9999984 3.248316e-06 1.624158e-06
[15,] 0.9999994 1.104936e-06 5.524681e-07
[16,] 1.0000000 2.237241e-08 1.118620e-08
[17,] 1.0000000 3.823168e-08 1.911584e-08
[18,] 1.0000000 4.073508e-10 2.036754e-10
[19,] 1.0000000 1.015765e-09 5.078824e-10
[20,] 1.0000000 1.490747e-10 7.453735e-11
[21,] 1.0000000 3.948052e-10 1.974026e-10
[22,] 1.0000000 4.444512e-10 2.222256e-10
[23,] 1.0000000 2.768040e-10 1.384020e-10
[24,] 1.0000000 7.147707e-10 3.573853e-10
[25,] 1.0000000 6.905850e-10 3.452925e-10
[26,] 1.0000000 4.028974e-10 2.014487e-10
[27,] 1.0000000 3.385988e-10 1.692994e-10
[28,] 1.0000000 5.707004e-10 2.853502e-10
[29,] 1.0000000 4.183891e-10 2.091946e-10
[30,] 1.0000000 7.324472e-10 3.662236e-10
[31,] 1.0000000 8.320268e-10 4.160134e-10
[32,] 1.0000000 2.057341e-09 1.028671e-09
[33,] 1.0000000 4.538981e-09 2.269491e-09
[34,] 1.0000000 2.215938e-09 1.107969e-09
[35,] 1.0000000 1.676306e-10 8.381528e-11
[36,] 1.0000000 2.092142e-12 1.046071e-12
[37,] 1.0000000 5.385114e-12 2.692557e-12
[38,] 1.0000000 1.397479e-11 6.987397e-12
[39,] 1.0000000 2.206494e-11 1.103247e-11
[40,] 1.0000000 3.089009e-11 1.544504e-11
[41,] 1.0000000 3.201753e-11 1.600877e-11
[42,] 1.0000000 3.952888e-11 1.976444e-11
[43,] 1.0000000 7.757441e-11 3.878721e-11
[44,] 1.0000000 2.070058e-10 1.035029e-10
[45,] 1.0000000 2.456057e-10 1.228029e-10
[46,] 1.0000000 3.306720e-11 1.653360e-11
[47,] 1.0000000 8.594485e-11 4.297242e-11
[48,] 1.0000000 9.768460e-11 4.884230e-11
[49,] 1.0000000 5.361826e-11 2.680913e-11
[50,] 1.0000000 7.445674e-11 3.722837e-11
[51,] 1.0000000 1.698065e-11 8.490324e-12
[52,] 1.0000000 1.313074e-11 6.565368e-12
[53,] 1.0000000 3.255383e-11 1.627692e-11
[54,] 1.0000000 2.928400e-11 1.464200e-11
[55,] 1.0000000 3.761371e-11 1.880686e-11
[56,] 1.0000000 1.188709e-10 5.943543e-11
[57,] 1.0000000 2.842452e-10 1.421226e-10
[58,] 1.0000000 4.063593e-10 2.031796e-10
[59,] 1.0000000 1.305428e-09 6.527141e-10
[60,] 1.0000000 6.646889e-11 3.323444e-11
[61,] 1.0000000 1.874980e-10 9.374899e-11
[62,] 1.0000000 5.014853e-10 2.507426e-10
[63,] 1.0000000 2.157230e-10 1.078615e-10
[64,] 1.0000000 7.728921e-11 3.864460e-11
[65,] 1.0000000 2.507362e-10 1.253681e-10
[66,] 1.0000000 9.569850e-11 4.784925e-11
[67,] 1.0000000 3.839132e-10 1.919566e-10
[68,] 1.0000000 9.864776e-10 4.932388e-10
[69,] 1.0000000 4.029839e-09 2.014920e-09
[70,] 1.0000000 1.499045e-08 7.495226e-09
[71,] 1.0000000 5.584652e-08 2.792326e-08
[72,] 0.9999999 1.220530e-07 6.102649e-08
[73,] 0.9999998 3.318007e-07 1.659004e-07
[74,] 0.9999996 7.685796e-07 3.842898e-07
[75,] 0.9999988 2.393849e-06 1.196924e-06
[76,] 0.9999974 5.193268e-06 2.596634e-06
[77,] 0.9999950 1.004255e-05 5.021274e-06
[78,] 0.9999967 6.558805e-06 3.279403e-06
[79,] 0.9999868 2.645540e-05 1.322770e-05
[80,] 0.9999474 1.051035e-04 5.255177e-05
[81,] 0.9999957 8.556326e-06 4.278163e-06
[82,] 0.9999787 4.253140e-05 2.126570e-05
[83,] 0.9999442 1.116141e-04 5.580704e-05
[84,] 0.9997077 5.846585e-04 2.923293e-04
[85,] 0.9985139 2.972259e-03 1.486130e-03
[86,] 0.9964060 7.188069e-03 3.594034e-03
[87,] 0.9993082 1.383571e-03 6.917854e-04
> postscript(file="/var/www/html/freestat/rcomp/tmp/10k2o1291209396.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/20k2o1291209396.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/30k2o1291209396.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/4ab1r1291209396.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/5ab1r1291209396.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 = 100
Frequency = 1
1 2 3 4 5 6
-32.3975016 246.4042470 -64.1825875 -32.2668349 -1.9126527 -88.9191456
7 8 9 10 11 12
23.1706089 1.0488270 -9.1375442 -4.1702078 -56.9617626 16.4778703
13 14 15 16 17 18
249.1643041 -66.8082221 107.2890111 -8.1938490 34.1640783 107.0810302
19 20 21 22 23 24
-3.6102466 149.6950004 -42.6155129 -94.2948308 -0.4973922 170.3706387
25 26 27 28 29 30
8.2573513 -64.1652179 48.4757564 -29.5666723 -52.1686897 25.0862936
31 32 33 34 35 36
-42.4844843 -38.8655831 -46.4078708 -31.4340446 -31.3255209 79.1756293
37 38 39 40 41 42
-18.2620557 3.0370182 24.0242092 68.8758297 120.7368368 -52.4553579
43 44 45 46 47 48
25.2690271 13.8808036 -26.3200290 -44.4004150 -21.4685826 62.6161384
49 50 51 52 53 54
-6.3259475 -12.5674051 23.4576746 65.8041459 13.8979125 8.0202868
55 56 57 58 59 60
42.2893766 22.2812051 -55.8567559 -21.1259901 3.3230343 -12.5344012
61 62 63 64 65 66
-35.8787864 -18.5885400 -41.9486996 -12.1405312 -2.2409532 144.5650848
67 68 69 70 71 72
-31.5257035 -17.4808560 26.0391452 13.8398924 -20.5414473 -31.0426578
73 74 75 76 77 78
-13.0630950 -36.0019600 -16.1877789 -41.0810988 -16.6718614 -35.4829494
79 80 81 82 83 84
3.0239725 -40.9464818 -20.8377252 -8.7335887 -10.9416332 -33.0183626
85 86 87 88 89 90
-28.7736231 -12.8257589 40.4195902 -30.8460046 -6.3334711 -30.0048113
91 92 93 94 95 96
-15.3429793 -26.6235452 14.0191418 -53.1653763 -10.3543609 -38.1392513
97 98 99 100
-33.1405226 -37.6969381 -33.5644456 -50.4118595
> postscript(file="/var/www/html/freestat/rcomp/tmp/6ab1r1291209396.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 = 100
Frequency = 1
lag(myerror, k = 1) myerror
0 -32.3975016 NA
1 246.4042470 -32.3975016
2 -64.1825875 246.4042470
3 -32.2668349 -64.1825875
4 -1.9126527 -32.2668349
5 -88.9191456 -1.9126527
6 23.1706089 -88.9191456
7 1.0488270 23.1706089
8 -9.1375442 1.0488270
9 -4.1702078 -9.1375442
10 -56.9617626 -4.1702078
11 16.4778703 -56.9617626
12 249.1643041 16.4778703
13 -66.8082221 249.1643041
14 107.2890111 -66.8082221
15 -8.1938490 107.2890111
16 34.1640783 -8.1938490
17 107.0810302 34.1640783
18 -3.6102466 107.0810302
19 149.6950004 -3.6102466
20 -42.6155129 149.6950004
21 -94.2948308 -42.6155129
22 -0.4973922 -94.2948308
23 170.3706387 -0.4973922
24 8.2573513 170.3706387
25 -64.1652179 8.2573513
26 48.4757564 -64.1652179
27 -29.5666723 48.4757564
28 -52.1686897 -29.5666723
29 25.0862936 -52.1686897
30 -42.4844843 25.0862936
31 -38.8655831 -42.4844843
32 -46.4078708 -38.8655831
33 -31.4340446 -46.4078708
34 -31.3255209 -31.4340446
35 79.1756293 -31.3255209
36 -18.2620557 79.1756293
37 3.0370182 -18.2620557
38 24.0242092 3.0370182
39 68.8758297 24.0242092
40 120.7368368 68.8758297
41 -52.4553579 120.7368368
42 25.2690271 -52.4553579
43 13.8808036 25.2690271
44 -26.3200290 13.8808036
45 -44.4004150 -26.3200290
46 -21.4685826 -44.4004150
47 62.6161384 -21.4685826
48 -6.3259475 62.6161384
49 -12.5674051 -6.3259475
50 23.4576746 -12.5674051
51 65.8041459 23.4576746
52 13.8979125 65.8041459
53 8.0202868 13.8979125
54 42.2893766 8.0202868
55 22.2812051 42.2893766
56 -55.8567559 22.2812051
57 -21.1259901 -55.8567559
58 3.3230343 -21.1259901
59 -12.5344012 3.3230343
60 -35.8787864 -12.5344012
61 -18.5885400 -35.8787864
62 -41.9486996 -18.5885400
63 -12.1405312 -41.9486996
64 -2.2409532 -12.1405312
65 144.5650848 -2.2409532
66 -31.5257035 144.5650848
67 -17.4808560 -31.5257035
68 26.0391452 -17.4808560
69 13.8398924 26.0391452
70 -20.5414473 13.8398924
71 -31.0426578 -20.5414473
72 -13.0630950 -31.0426578
73 -36.0019600 -13.0630950
74 -16.1877789 -36.0019600
75 -41.0810988 -16.1877789
76 -16.6718614 -41.0810988
77 -35.4829494 -16.6718614
78 3.0239725 -35.4829494
79 -40.9464818 3.0239725
80 -20.8377252 -40.9464818
81 -8.7335887 -20.8377252
82 -10.9416332 -8.7335887
83 -33.0183626 -10.9416332
84 -28.7736231 -33.0183626
85 -12.8257589 -28.7736231
86 40.4195902 -12.8257589
87 -30.8460046 40.4195902
88 -6.3334711 -30.8460046
89 -30.0048113 -6.3334711
90 -15.3429793 -30.0048113
91 -26.6235452 -15.3429793
92 14.0191418 -26.6235452
93 -53.1653763 14.0191418
94 -10.3543609 -53.1653763
95 -38.1392513 -10.3543609
96 -33.1405226 -38.1392513
97 -37.6969381 -33.1405226
98 -33.5644456 -37.6969381
99 -50.4118595 -33.5644456
100 NA -50.4118595
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 246.4042470 -32.3975016
[2,] -64.1825875 246.4042470
[3,] -32.2668349 -64.1825875
[4,] -1.9126527 -32.2668349
[5,] -88.9191456 -1.9126527
[6,] 23.1706089 -88.9191456
[7,] 1.0488270 23.1706089
[8,] -9.1375442 1.0488270
[9,] -4.1702078 -9.1375442
[10,] -56.9617626 -4.1702078
[11,] 16.4778703 -56.9617626
[12,] 249.1643041 16.4778703
[13,] -66.8082221 249.1643041
[14,] 107.2890111 -66.8082221
[15,] -8.1938490 107.2890111
[16,] 34.1640783 -8.1938490
[17,] 107.0810302 34.1640783
[18,] -3.6102466 107.0810302
[19,] 149.6950004 -3.6102466
[20,] -42.6155129 149.6950004
[21,] -94.2948308 -42.6155129
[22,] -0.4973922 -94.2948308
[23,] 170.3706387 -0.4973922
[24,] 8.2573513 170.3706387
[25,] -64.1652179 8.2573513
[26,] 48.4757564 -64.1652179
[27,] -29.5666723 48.4757564
[28,] -52.1686897 -29.5666723
[29,] 25.0862936 -52.1686897
[30,] -42.4844843 25.0862936
[31,] -38.8655831 -42.4844843
[32,] -46.4078708 -38.8655831
[33,] -31.4340446 -46.4078708
[34,] -31.3255209 -31.4340446
[35,] 79.1756293 -31.3255209
[36,] -18.2620557 79.1756293
[37,] 3.0370182 -18.2620557
[38,] 24.0242092 3.0370182
[39,] 68.8758297 24.0242092
[40,] 120.7368368 68.8758297
[41,] -52.4553579 120.7368368
[42,] 25.2690271 -52.4553579
[43,] 13.8808036 25.2690271
[44,] -26.3200290 13.8808036
[45,] -44.4004150 -26.3200290
[46,] -21.4685826 -44.4004150
[47,] 62.6161384 -21.4685826
[48,] -6.3259475 62.6161384
[49,] -12.5674051 -6.3259475
[50,] 23.4576746 -12.5674051
[51,] 65.8041459 23.4576746
[52,] 13.8979125 65.8041459
[53,] 8.0202868 13.8979125
[54,] 42.2893766 8.0202868
[55,] 22.2812051 42.2893766
[56,] -55.8567559 22.2812051
[57,] -21.1259901 -55.8567559
[58,] 3.3230343 -21.1259901
[59,] -12.5344012 3.3230343
[60,] -35.8787864 -12.5344012
[61,] -18.5885400 -35.8787864
[62,] -41.9486996 -18.5885400
[63,] -12.1405312 -41.9486996
[64,] -2.2409532 -12.1405312
[65,] 144.5650848 -2.2409532
[66,] -31.5257035 144.5650848
[67,] -17.4808560 -31.5257035
[68,] 26.0391452 -17.4808560
[69,] 13.8398924 26.0391452
[70,] -20.5414473 13.8398924
[71,] -31.0426578 -20.5414473
[72,] -13.0630950 -31.0426578
[73,] -36.0019600 -13.0630950
[74,] -16.1877789 -36.0019600
[75,] -41.0810988 -16.1877789
[76,] -16.6718614 -41.0810988
[77,] -35.4829494 -16.6718614
[78,] 3.0239725 -35.4829494
[79,] -40.9464818 3.0239725
[80,] -20.8377252 -40.9464818
[81,] -8.7335887 -20.8377252
[82,] -10.9416332 -8.7335887
[83,] -33.0183626 -10.9416332
[84,] -28.7736231 -33.0183626
[85,] -12.8257589 -28.7736231
[86,] 40.4195902 -12.8257589
[87,] -30.8460046 40.4195902
[88,] -6.3334711 -30.8460046
[89,] -30.0048113 -6.3334711
[90,] -15.3429793 -30.0048113
[91,] -26.6235452 -15.3429793
[92,] 14.0191418 -26.6235452
[93,] -53.1653763 14.0191418
[94,] -10.3543609 -53.1653763
[95,] -38.1392513 -10.3543609
[96,] -33.1405226 -38.1392513
[97,] -37.6969381 -33.1405226
[98,] -33.5644456 -37.6969381
[99,] -50.4118595 -33.5644456
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 246.4042470 -32.3975016
2 -64.1825875 246.4042470
3 -32.2668349 -64.1825875
4 -1.9126527 -32.2668349
5 -88.9191456 -1.9126527
6 23.1706089 -88.9191456
7 1.0488270 23.1706089
8 -9.1375442 1.0488270
9 -4.1702078 -9.1375442
10 -56.9617626 -4.1702078
11 16.4778703 -56.9617626
12 249.1643041 16.4778703
13 -66.8082221 249.1643041
14 107.2890111 -66.8082221
15 -8.1938490 107.2890111
16 34.1640783 -8.1938490
17 107.0810302 34.1640783
18 -3.6102466 107.0810302
19 149.6950004 -3.6102466
20 -42.6155129 149.6950004
21 -94.2948308 -42.6155129
22 -0.4973922 -94.2948308
23 170.3706387 -0.4973922
24 8.2573513 170.3706387
25 -64.1652179 8.2573513
26 48.4757564 -64.1652179
27 -29.5666723 48.4757564
28 -52.1686897 -29.5666723
29 25.0862936 -52.1686897
30 -42.4844843 25.0862936
31 -38.8655831 -42.4844843
32 -46.4078708 -38.8655831
33 -31.4340446 -46.4078708
34 -31.3255209 -31.4340446
35 79.1756293 -31.3255209
36 -18.2620557 79.1756293
37 3.0370182 -18.2620557
38 24.0242092 3.0370182
39 68.8758297 24.0242092
40 120.7368368 68.8758297
41 -52.4553579 120.7368368
42 25.2690271 -52.4553579
43 13.8808036 25.2690271
44 -26.3200290 13.8808036
45 -44.4004150 -26.3200290
46 -21.4685826 -44.4004150
47 62.6161384 -21.4685826
48 -6.3259475 62.6161384
49 -12.5674051 -6.3259475
50 23.4576746 -12.5674051
51 65.8041459 23.4576746
52 13.8979125 65.8041459
53 8.0202868 13.8979125
54 42.2893766 8.0202868
55 22.2812051 42.2893766
56 -55.8567559 22.2812051
57 -21.1259901 -55.8567559
58 3.3230343 -21.1259901
59 -12.5344012 3.3230343
60 -35.8787864 -12.5344012
61 -18.5885400 -35.8787864
62 -41.9486996 -18.5885400
63 -12.1405312 -41.9486996
64 -2.2409532 -12.1405312
65 144.5650848 -2.2409532
66 -31.5257035 144.5650848
67 -17.4808560 -31.5257035
68 26.0391452 -17.4808560
69 13.8398924 26.0391452
70 -20.5414473 13.8398924
71 -31.0426578 -20.5414473
72 -13.0630950 -31.0426578
73 -36.0019600 -13.0630950
74 -16.1877789 -36.0019600
75 -41.0810988 -16.1877789
76 -16.6718614 -41.0810988
77 -35.4829494 -16.6718614
78 3.0239725 -35.4829494
79 -40.9464818 3.0239725
80 -20.8377252 -40.9464818
81 -8.7335887 -20.8377252
82 -10.9416332 -8.7335887
83 -33.0183626 -10.9416332
84 -28.7736231 -33.0183626
85 -12.8257589 -28.7736231
86 40.4195902 -12.8257589
87 -30.8460046 40.4195902
88 -6.3334711 -30.8460046
89 -30.0048113 -6.3334711
90 -15.3429793 -30.0048113
91 -26.6235452 -15.3429793
92 14.0191418 -26.6235452
93 -53.1653763 14.0191418
94 -10.3543609 -53.1653763
95 -38.1392513 -10.3543609
96 -33.1405226 -38.1392513
97 -37.6969381 -33.1405226
98 -33.5644456 -37.6969381
99 -50.4118595 -33.5644456
> 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/7lk1u1291209396.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/8ebif1291209396.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/freestat/rcomp/tmp/9ebif1291209396.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> opar <- par(mfrow = c(2,2), oma = c(0, 0, 1.1, 0))
> plot(mylm, las = 1, sub='Residual Diagnostics')
> par(opar)
> dev.off()
null device
1
> if (n > n25) {
+ postscript(file="/var/www/html/freestat/rcomp/tmp/10ebif1291209396.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/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/11a3go1291209396.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/12luxr1291209396.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/13rdc31291209396.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/14kntn1291209396.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/1555st1291209396.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/162x721291209396.tab")
+ }
>
> try(system("convert tmp/10k2o1291209396.ps tmp/10k2o1291209396.png",intern=TRUE))
character(0)
> try(system("convert tmp/20k2o1291209396.ps tmp/20k2o1291209396.png",intern=TRUE))
character(0)
> try(system("convert tmp/30k2o1291209396.ps tmp/30k2o1291209396.png",intern=TRUE))
character(0)
> try(system("convert tmp/4ab1r1291209396.ps tmp/4ab1r1291209396.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ab1r1291209396.ps tmp/5ab1r1291209396.png",intern=TRUE))
character(0)
> try(system("convert tmp/6ab1r1291209396.ps tmp/6ab1r1291209396.png",intern=TRUE))
character(0)
> try(system("convert tmp/7lk1u1291209396.ps tmp/7lk1u1291209396.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ebif1291209396.ps tmp/8ebif1291209396.png",intern=TRUE))
character(0)
> try(system("convert tmp/9ebif1291209396.ps tmp/9ebif1291209396.png",intern=TRUE))
character(0)
> try(system("convert tmp/10ebif1291209396.ps tmp/10ebif1291209396.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
4.454 2.567 4.769