R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> x <- array(list(16203
+ ,112
+ ,13808
+ ,11752
+ ,10751
+ ,10144
+ ,17432
+ ,112
+ ,16203
+ ,13808
+ ,11752
+ ,10751
+ ,18014
+ ,304
+ ,17432
+ ,16203
+ ,13808
+ ,11752
+ ,16956
+ ,794
+ ,18014
+ ,17432
+ ,16203
+ ,13808
+ ,17982
+ ,901
+ ,16956
+ ,18014
+ ,17432
+ ,16203
+ ,19435
+ ,1232
+ ,17982
+ ,16956
+ ,18014
+ ,17432
+ ,19990
+ ,1240
+ ,19435
+ ,17982
+ ,16956
+ ,18014
+ ,20154
+ ,1032
+ ,19990
+ ,19435
+ ,17982
+ ,16956
+ ,10327
+ ,1145
+ ,20154
+ ,19990
+ ,19435
+ ,17982
+ ,9807
+ ,1588
+ ,10327
+ ,20154
+ ,19990
+ ,19435
+ ,10862
+ ,2264
+ ,9807
+ ,10327
+ ,20154
+ ,19990
+ ,13743
+ ,2209
+ ,10862
+ ,9807
+ ,10327
+ ,20154
+ ,16458
+ ,2917
+ ,13743
+ ,10862
+ ,9807
+ ,10144
+ ,18466
+ ,243
+ ,16458
+ ,13743
+ ,10862
+ ,10751
+ ,18810
+ ,558
+ ,18466
+ ,16458
+ ,13743
+ ,11752
+ ,17361
+ ,1238
+ ,18810
+ ,18466
+ ,16458
+ ,13808
+ ,17411
+ ,1502
+ ,17361
+ ,18810
+ ,18466
+ ,16203
+ ,18517
+ ,2000
+ ,17411
+ ,17361
+ ,18810
+ ,17432
+ ,18525
+ ,2146
+ ,18517
+ ,17411
+ ,17361
+ ,18014
+ ,17859
+ ,2066
+ ,18525
+ ,18517
+ ,17411
+ ,16956
+ ,9499
+ ,2046
+ ,17859
+ ,18525
+ ,18517
+ ,17982
+ ,9490
+ ,1952
+ ,9499
+ ,17859
+ ,18525
+ ,19435
+ ,9255
+ ,2771
+ ,9490
+ ,9499
+ ,17859
+ ,19990
+ ,10758
+ ,3278
+ ,9255
+ ,9490
+ ,9499
+ ,20154
+ ,12375
+ ,4000
+ ,10758
+ ,9255
+ ,9490
+ ,10327
+ ,14617
+ ,410
+ ,12375
+ ,10758
+ ,9255
+ ,9807
+ ,15427
+ ,1107
+ ,14617
+ ,12375
+ ,10758
+ ,10862
+ ,14136
+ ,1622
+ ,15427
+ ,14617
+ ,12375
+ ,13743
+ ,14308
+ ,1986
+ ,14136
+ ,15427
+ ,14617
+ ,16458
+ ,15293
+ ,2036
+ ,14308
+ ,14136
+ ,15427
+ ,18466
+ ,15679
+ ,2400
+ ,15293
+ ,14308
+ ,14136
+ ,18810
+ ,16319
+ ,2736
+ ,15679
+ ,15293
+ ,14308
+ ,17361
+ ,11196
+ ,2901
+ ,16319
+ ,15679
+ ,15293
+ ,17411
+ ,11169
+ ,2883
+ ,11196
+ ,16319
+ ,15679
+ ,18517
+ ,12158
+ ,3747
+ ,11169
+ ,11196
+ ,16319
+ ,18525
+ ,14251
+ ,4075
+ ,12158
+ ,11169
+ ,11196
+ ,17859
+ ,16237
+ ,4996
+ ,14251
+ ,12158
+ ,11169
+ ,9499
+ ,19706
+ ,575
+ ,16237
+ ,14251
+ ,12158
+ ,9490
+ ,18960
+ ,999
+ ,19706
+ ,16237
+ ,14251
+ ,9255
+ ,18537
+ ,1411
+ ,18960
+ ,19706
+ ,16237
+ ,10758
+ ,19103
+ ,1493
+ ,18537
+ ,18960
+ ,19706
+ ,12375
+ ,19691
+ ,1846
+ ,19103
+ ,18537
+ ,18960
+ ,14617
+ ,19464
+ ,2899
+ ,19691
+ ,19103
+ ,18537
+ ,15427
+ ,17264
+ ,2372
+ ,19464
+ ,19691
+ ,19103
+ ,14136
+ ,8957
+ ,2856
+ ,17264
+ ,19464
+ ,19691
+ ,14308
+ ,9703
+ ,3468
+ ,8957
+ ,17264
+ ,19464
+ ,15293
+ ,9166
+ ,4193
+ ,9703
+ ,8957
+ ,17264
+ ,15679
+ ,9519
+ ,4440
+ ,9166
+ ,9703
+ ,8957
+ ,16319
+ ,10535
+ ,4186
+ ,9519
+ ,9166
+ ,9703
+ ,11196
+ ,11526
+ ,655
+ ,10535
+ ,9519
+ ,9166
+ ,11169
+ ,9630
+ ,1453
+ ,11526
+ ,10535
+ ,9519
+ ,12158
+ ,7061
+ ,1989
+ ,9630
+ ,11526
+ ,10535
+ ,14251
+ ,6021
+ ,2209
+ ,7061
+ ,9630
+ ,11526
+ ,16237)
+ ,dim=c(6
+ ,53)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y1'
+ ,'Y2'
+ ,'Y3'
+ ,'Y4')
+ ,1:53))
> y <- array(NA,dim=c(6,53),dimnames=list(c('Y','X','Y1','Y2','Y3','Y4'),1:53))
> 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 Y1 Y2 Y3 Y4 M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 t
1 16203 112 13808 11752 10751 10144 1 0 0 0 0 0 0 0 0 0 0 1
2 17432 112 16203 13808 11752 10751 0 1 0 0 0 0 0 0 0 0 0 2
3 18014 304 17432 16203 13808 11752 0 0 1 0 0 0 0 0 0 0 0 3
4 16956 794 18014 17432 16203 13808 0 0 0 1 0 0 0 0 0 0 0 4
5 17982 901 16956 18014 17432 16203 0 0 0 0 1 0 0 0 0 0 0 5
6 19435 1232 17982 16956 18014 17432 0 0 0 0 0 1 0 0 0 0 0 6
7 19990 1240 19435 17982 16956 18014 0 0 0 0 0 0 1 0 0 0 0 7
8 20154 1032 19990 19435 17982 16956 0 0 0 0 0 0 0 1 0 0 0 8
9 10327 1145 20154 19990 19435 17982 0 0 0 0 0 0 0 0 1 0 0 9
10 9807 1588 10327 20154 19990 19435 0 0 0 0 0 0 0 0 0 1 0 10
11 10862 2264 9807 10327 20154 19990 0 0 0 0 0 0 0 0 0 0 1 11
12 13743 2209 10862 9807 10327 20154 0 0 0 0 0 0 0 0 0 0 0 12
13 16458 2917 13743 10862 9807 10144 1 0 0 0 0 0 0 0 0 0 0 13
14 18466 243 16458 13743 10862 10751 0 1 0 0 0 0 0 0 0 0 0 14
15 18810 558 18466 16458 13743 11752 0 0 1 0 0 0 0 0 0 0 0 15
16 17361 1238 18810 18466 16458 13808 0 0 0 1 0 0 0 0 0 0 0 16
17 17411 1502 17361 18810 18466 16203 0 0 0 0 1 0 0 0 0 0 0 17
18 18517 2000 17411 17361 18810 17432 0 0 0 0 0 1 0 0 0 0 0 18
19 18525 2146 18517 17411 17361 18014 0 0 0 0 0 0 1 0 0 0 0 19
20 17859 2066 18525 18517 17411 16956 0 0 0 0 0 0 0 1 0 0 0 20
21 9499 2046 17859 18525 18517 17982 0 0 0 0 0 0 0 0 1 0 0 21
22 9490 1952 9499 17859 18525 19435 0 0 0 0 0 0 0 0 0 1 0 22
23 9255 2771 9490 9499 17859 19990 0 0 0 0 0 0 0 0 0 0 1 23
24 10758 3278 9255 9490 9499 20154 0 0 0 0 0 0 0 0 0 0 0 24
25 12375 4000 10758 9255 9490 10327 1 0 0 0 0 0 0 0 0 0 0 25
26 14617 410 12375 10758 9255 9807 0 1 0 0 0 0 0 0 0 0 0 26
27 15427 1107 14617 12375 10758 10862 0 0 1 0 0 0 0 0 0 0 0 27
28 14136 1622 15427 14617 12375 13743 0 0 0 1 0 0 0 0 0 0 0 28
29 14308 1986 14136 15427 14617 16458 0 0 0 0 1 0 0 0 0 0 0 29
30 15293 2036 14308 14136 15427 18466 0 0 0 0 0 1 0 0 0 0 0 30
31 15679 2400 15293 14308 14136 18810 0 0 0 0 0 0 1 0 0 0 0 31
32 16319 2736 15679 15293 14308 17361 0 0 0 0 0 0 0 1 0 0 0 32
33 11196 2901 16319 15679 15293 17411 0 0 0 0 0 0 0 0 1 0 0 33
34 11169 2883 11196 16319 15679 18517 0 0 0 0 0 0 0 0 0 1 0 34
35 12158 3747 11169 11196 16319 18525 0 0 0 0 0 0 0 0 0 0 1 35
36 14251 4075 12158 11169 11196 17859 0 0 0 0 0 0 0 0 0 0 0 36
37 16237 4996 14251 12158 11169 9499 1 0 0 0 0 0 0 0 0 0 0 37
38 19706 575 16237 14251 12158 9490 0 1 0 0 0 0 0 0 0 0 0 38
39 18960 999 19706 16237 14251 9255 0 0 1 0 0 0 0 0 0 0 0 39
40 18537 1411 18960 19706 16237 10758 0 0 0 1 0 0 0 0 0 0 0 40
41 19103 1493 18537 18960 19706 12375 0 0 0 0 1 0 0 0 0 0 0 41
42 19691 1846 19103 18537 18960 14617 0 0 0 0 0 1 0 0 0 0 0 42
43 19464 2899 19691 19103 18537 15427 0 0 0 0 0 0 1 0 0 0 0 43
44 17264 2372 19464 19691 19103 14136 0 0 0 0 0 0 0 1 0 0 0 44
45 8957 2856 17264 19464 19691 14308 0 0 0 0 0 0 0 0 1 0 0 45
46 9703 3468 8957 17264 19464 15293 0 0 0 0 0 0 0 0 0 1 0 46
47 9166 4193 9703 8957 17264 15679 0 0 0 0 0 0 0 0 0 0 1 47
48 9519 4440 9166 9703 8957 16319 0 0 0 0 0 0 0 0 0 0 0 48
49 10535 4186 9519 9166 9703 11196 1 0 0 0 0 0 0 0 0 0 0 49
50 11526 655 10535 9519 9166 11169 0 1 0 0 0 0 0 0 0 0 0 50
51 9630 1453 11526 10535 9519 12158 0 0 1 0 0 0 0 0 0 0 0 51
52 7061 1989 9630 11526 10535 14251 0 0 0 1 0 0 0 0 0 0 0 52
53 6021 2209 7061 9630 11526 16237 0 0 0 0 1 0 0 0 0 0 0 53
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X Y1 Y2 Y3 Y4
5039.5227 0.1526 1.2172 -0.1589 -0.1334 -0.1251
M1 M2 M3 M4 M5 M6
-1237.4867 -796.2316 -2800.8395 -3149.6993 -830.7427 -166.8952
M7 M8 M9 M10 M11 t
-1271.2732 -1886.8124 -8925.7116 834.0082 -272.6532 -28.6099
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2217.64 -347.00 -51.34 348.39 2429.17
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5039.5227 4459.3614 1.130 0.26612
X 0.1526 0.2696 0.566 0.57510
Y1 1.2172 0.1774 6.861 5.78e-08 ***
Y2 -0.1589 0.2690 -0.591 0.55847
Y3 -0.1334 0.1932 -0.690 0.49458
Y4 -0.1251 0.1681 -0.744 0.46166
M1 -1237.4867 1506.3151 -0.822 0.41690
M2 -796.2316 1599.2779 -0.498 0.62169
M3 -2800.8395 1352.0242 -2.072 0.04574 *
M4 -3149.6993 1070.0869 -2.943 0.00573 **
M5 -830.7427 944.2311 -0.880 0.38496
M6 -166.8952 1010.6006 -0.165 0.86978
M7 -1271.2732 970.7240 -1.310 0.19886
M8 -1886.8124 952.8290 -1.980 0.05558 .
M9 -8925.7116 1009.2617 -8.844 1.91e-10 ***
M10 834.0082 1625.6221 0.513 0.61115
M11 -272.6532 1649.8824 -0.165 0.86969
t -28.6099 19.3111 -1.482 0.14741
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 937.1 on 35 degrees of freedom
Multiple R-squared: 0.9634, Adjusted R-squared: 0.9456
F-statistic: 54.16 on 17 and 35 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.26037246 0.5207449 0.7396275
[2,] 0.12558574 0.2511715 0.8744143
[3,] 0.28800041 0.5760008 0.7119996
[4,] 0.21497285 0.4299457 0.7850272
[5,] 0.12533178 0.2506636 0.8746682
[6,] 0.09968680 0.1993736 0.9003132
[7,] 0.06053384 0.1210677 0.9394662
[8,] 0.05325616 0.1065123 0.9467438
[9,] 0.05116029 0.1023206 0.9488397
[10,] 0.05245655 0.1049131 0.9475434
[11,] 0.17485183 0.3497037 0.8251482
[12,] 0.22338971 0.4467794 0.7766103
> postscript(file="/var/www/html/rcomp/tmp/163ay1258725288.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(x[,1], type='l', main='Actuals and Interpolation', ylab='value of Actuals and Interpolation (dots)', xlab='time or index')
> points(x[,1]-mysum$resid)
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/2nb321258725288.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> plot(mysum$resid, type='b', pch=19, main='Residuals', ylab='value of Residuals', xlab='time or index')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/3b7l71258725288.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> hist(mysum$resid, main='Residual Histogram', xlab='values of Residuals')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/43ox61258725288.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> densityplot(~mysum$resid,col='black',main='Residual Density Plot', xlab='values of Residuals')
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/5e5jx1258725288.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 = 53
Frequency = 1
1 2 3 4 5 6
175.525980 -1387.184238 482.809881 -208.987865 354.137351 -64.239814
7 8 9 10 11 12
-51.343964 348.386790 -2217.641745 -292.947497 956.778349 945.200859
13 14 15 16 17 18
157.670926 -469.271676 356.631483 -298.983923 -193.806686 109.459445
19 20 21 22 23 24
-230.575299 -199.871749 -401.527737 125.590044 -436.068058 -64.321976
25 26 27 28 29 30
-388.623889 162.680349 759.954224 -285.768422 -120.876808 165.942366
31 32 33 34 35 36
328.620168 1089.882109 2429.167206 -798.878420 498.676057 322.904339
37 38 39 40 41 42
-5.301635 1771.561498 -663.067289 1140.742252 465.147365 -211.161996
43 44 45 46 47 48
-46.700905 -1238.397150 190.002275 966.235874 -1019.386348 -1203.783222
49 50 51 52 53
60.728617 -77.785933 -936.328298 -347.002041 -504.601221
> postscript(file="/var/www/html/rcomp/tmp/6yusr1258725288.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 = 53
Frequency = 1
lag(myerror, k = 1) myerror
0 175.525980 NA
1 -1387.184238 175.525980
2 482.809881 -1387.184238
3 -208.987865 482.809881
4 354.137351 -208.987865
5 -64.239814 354.137351
6 -51.343964 -64.239814
7 348.386790 -51.343964
8 -2217.641745 348.386790
9 -292.947497 -2217.641745
10 956.778349 -292.947497
11 945.200859 956.778349
12 157.670926 945.200859
13 -469.271676 157.670926
14 356.631483 -469.271676
15 -298.983923 356.631483
16 -193.806686 -298.983923
17 109.459445 -193.806686
18 -230.575299 109.459445
19 -199.871749 -230.575299
20 -401.527737 -199.871749
21 125.590044 -401.527737
22 -436.068058 125.590044
23 -64.321976 -436.068058
24 -388.623889 -64.321976
25 162.680349 -388.623889
26 759.954224 162.680349
27 -285.768422 759.954224
28 -120.876808 -285.768422
29 165.942366 -120.876808
30 328.620168 165.942366
31 1089.882109 328.620168
32 2429.167206 1089.882109
33 -798.878420 2429.167206
34 498.676057 -798.878420
35 322.904339 498.676057
36 -5.301635 322.904339
37 1771.561498 -5.301635
38 -663.067289 1771.561498
39 1140.742252 -663.067289
40 465.147365 1140.742252
41 -211.161996 465.147365
42 -46.700905 -211.161996
43 -1238.397150 -46.700905
44 190.002275 -1238.397150
45 966.235874 190.002275
46 -1019.386348 966.235874
47 -1203.783222 -1019.386348
48 60.728617 -1203.783222
49 -77.785933 60.728617
50 -936.328298 -77.785933
51 -347.002041 -936.328298
52 -504.601221 -347.002041
53 NA -504.601221
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] -1387.184238 175.525980
[2,] 482.809881 -1387.184238
[3,] -208.987865 482.809881
[4,] 354.137351 -208.987865
[5,] -64.239814 354.137351
[6,] -51.343964 -64.239814
[7,] 348.386790 -51.343964
[8,] -2217.641745 348.386790
[9,] -292.947497 -2217.641745
[10,] 956.778349 -292.947497
[11,] 945.200859 956.778349
[12,] 157.670926 945.200859
[13,] -469.271676 157.670926
[14,] 356.631483 -469.271676
[15,] -298.983923 356.631483
[16,] -193.806686 -298.983923
[17,] 109.459445 -193.806686
[18,] -230.575299 109.459445
[19,] -199.871749 -230.575299
[20,] -401.527737 -199.871749
[21,] 125.590044 -401.527737
[22,] -436.068058 125.590044
[23,] -64.321976 -436.068058
[24,] -388.623889 -64.321976
[25,] 162.680349 -388.623889
[26,] 759.954224 162.680349
[27,] -285.768422 759.954224
[28,] -120.876808 -285.768422
[29,] 165.942366 -120.876808
[30,] 328.620168 165.942366
[31,] 1089.882109 328.620168
[32,] 2429.167206 1089.882109
[33,] -798.878420 2429.167206
[34,] 498.676057 -798.878420
[35,] 322.904339 498.676057
[36,] -5.301635 322.904339
[37,] 1771.561498 -5.301635
[38,] -663.067289 1771.561498
[39,] 1140.742252 -663.067289
[40,] 465.147365 1140.742252
[41,] -211.161996 465.147365
[42,] -46.700905 -211.161996
[43,] -1238.397150 -46.700905
[44,] 190.002275 -1238.397150
[45,] 966.235874 190.002275
[46,] -1019.386348 966.235874
[47,] -1203.783222 -1019.386348
[48,] 60.728617 -1203.783222
[49,] -77.785933 60.728617
[50,] -936.328298 -77.785933
[51,] -347.002041 -936.328298
[52,] -504.601221 -347.002041
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 -1387.184238 175.525980
2 482.809881 -1387.184238
3 -208.987865 482.809881
4 354.137351 -208.987865
5 -64.239814 354.137351
6 -51.343964 -64.239814
7 348.386790 -51.343964
8 -2217.641745 348.386790
9 -292.947497 -2217.641745
10 956.778349 -292.947497
11 945.200859 956.778349
12 157.670926 945.200859
13 -469.271676 157.670926
14 356.631483 -469.271676
15 -298.983923 356.631483
16 -193.806686 -298.983923
17 109.459445 -193.806686
18 -230.575299 109.459445
19 -199.871749 -230.575299
20 -401.527737 -199.871749
21 125.590044 -401.527737
22 -436.068058 125.590044
23 -64.321976 -436.068058
24 -388.623889 -64.321976
25 162.680349 -388.623889
26 759.954224 162.680349
27 -285.768422 759.954224
28 -120.876808 -285.768422
29 165.942366 -120.876808
30 328.620168 165.942366
31 1089.882109 328.620168
32 2429.167206 1089.882109
33 -798.878420 2429.167206
34 498.676057 -798.878420
35 322.904339 498.676057
36 -5.301635 322.904339
37 1771.561498 -5.301635
38 -663.067289 1771.561498
39 1140.742252 -663.067289
40 465.147365 1140.742252
41 -211.161996 465.147365
42 -46.700905 -211.161996
43 -1238.397150 -46.700905
44 190.002275 -1238.397150
45 966.235874 190.002275
46 -1019.386348 966.235874
47 -1203.783222 -1019.386348
48 60.728617 -1203.783222
49 -77.785933 60.728617
50 -936.328298 -77.785933
51 -347.002041 -936.328298
52 -504.601221 -347.002041
> plot(z,main=paste('Residual Lag plot, lowess, and regression line'), ylab='values of Residuals', xlab='lagged values of Residuals')
> lines(lowess(z))
> abline(lm(z))
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/75my21258725288.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> acf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/8kcso1258725288.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
> pacf(mysum$resid, lag.max=length(mysum$resid)/2, main='Residual Partial Autocorrelation Function')
> grid()
> dev.off()
null device
1
> postscript(file="/var/www/html/rcomp/tmp/9ui2d1258725288.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/rcomp/tmp/10u5lr1258725288.ps",horizontal=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556)
+ plot(kp3:nmkm3,gqarr[,2], main='Goldfeld-Quandt test',ylab='2-sided p-value',xlab='breakpoint')
+ grid()
+ dev.off()
+ }
null device
1
>
> #Note: the /var/www/html/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/www/html/rcomp/createtable")
>
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Estimated Regression Equation', 1, TRUE)
> a<-table.row.end(a)
> myeq <- colnames(x)[1]
> myeq <- paste(myeq, '[t] = ', sep='')
> for (i in 1:k){
+ if (mysum$coefficients[i,1] > 0) myeq <- paste(myeq, '+', '')
+ myeq <- paste(myeq, mysum$coefficients[i,1], sep=' ')
+ if (rownames(mysum$coefficients)[i] != '(Intercept)') {
+ myeq <- paste(myeq, rownames(mysum$coefficients)[i], sep='')
+ if (rownames(mysum$coefficients)[i] != 't') myeq <- paste(myeq, '[t]', sep='')
+ }
+ }
> myeq <- paste(myeq, ' + e[t]')
> a<-table.row.start(a)
> a<-table.element(a, myeq)
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/11lpbw1258725288.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a,hyperlink('http://www.xycoon.com/ols1.htm','Multiple Linear Regression - Ordinary Least Squares',''), 6, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a,'Variable',header=TRUE)
> a<-table.element(a,'Parameter',header=TRUE)
> a<-table.element(a,'S.D.',header=TRUE)
> a<-table.element(a,'T-STAT
H0: parameter = 0',header=TRUE)
> a<-table.element(a,'2-tail p-value',header=TRUE)
> a<-table.element(a,'1-tail p-value',header=TRUE)
> a<-table.row.end(a)
> for (i in 1:k){
+ a<-table.row.start(a)
+ a<-table.element(a,rownames(mysum$coefficients)[i],header=TRUE)
+ a<-table.element(a,mysum$coefficients[i,1])
+ a<-table.element(a, round(mysum$coefficients[i,2],6))
+ a<-table.element(a, round(mysum$coefficients[i,3],4))
+ a<-table.element(a, round(mysum$coefficients[i,4],6))
+ a<-table.element(a, round(mysum$coefficients[i,4]/2,6))
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/12lowq1258725288.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Regression Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple R',1,TRUE)
> a<-table.element(a, sqrt(mysum$r.squared))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'R-squared',1,TRUE)
> a<-table.element(a, mysum$r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Adjusted R-squared',1,TRUE)
> a<-table.element(a, mysum$adj.r.squared)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (value)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[1])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF numerator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[2])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'F-TEST (DF denominator)',1,TRUE)
> a<-table.element(a, mysum$fstatistic[3])
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'p-value',1,TRUE)
> a<-table.element(a, 1-pf(mysum$fstatistic[1],mysum$fstatistic[2],mysum$fstatistic[3]))
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Residual Statistics', 2, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Residual Standard Deviation',1,TRUE)
> a<-table.element(a, mysum$sigma)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Sum Squared Residuals',1,TRUE)
> a<-table.element(a, sum(myerror*myerror))
> a<-table.row.end(a)
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/13r3e91258725288.tab")
> a<-table.start()
> a<-table.row.start(a)
> a<-table.element(a, 'Multiple Linear Regression - Actuals, Interpolation, and Residuals', 4, TRUE)
> a<-table.row.end(a)
> a<-table.row.start(a)
> a<-table.element(a, 'Time or Index', 1, TRUE)
> a<-table.element(a, 'Actuals', 1, TRUE)
> a<-table.element(a, 'Interpolation
Forecast', 1, TRUE)
> a<-table.element(a, 'Residuals
Prediction Error', 1, TRUE)
> a<-table.row.end(a)
> for (i in 1:n) {
+ a<-table.row.start(a)
+ a<-table.element(a,i, 1, TRUE)
+ a<-table.element(a,x[i])
+ a<-table.element(a,x[i]-mysum$resid[i])
+ a<-table.element(a,mysum$resid[i])
+ a<-table.row.end(a)
+ }
> a<-table.end(a)
> table.save(a,file="/var/www/html/rcomp/tmp/14anjw1258725288.tab")
> if (n > n25) {
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'p-values',header=TRUE)
+ a<-table.element(a,'Alternative Hypothesis',3,header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'breakpoint index',header=TRUE)
+ a<-table.element(a,'greater',header=TRUE)
+ a<-table.element(a,'2-sided',header=TRUE)
+ a<-table.element(a,'less',header=TRUE)
+ a<-table.row.end(a)
+ for (mypoint in kp3:nmkm3) {
+ a<-table.row.start(a)
+ a<-table.element(a,mypoint,header=TRUE)
+ a<-table.element(a,gqarr[mypoint-kp3+1,1])
+ a<-table.element(a,gqarr[mypoint-kp3+1,2])
+ a<-table.element(a,gqarr[mypoint-kp3+1,3])
+ a<-table.row.end(a)
+ }
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/15tqp91258725288.tab")
+ a<-table.start()
+ a<-table.row.start(a)
+ a<-table.element(a,'Meta Analysis of Goldfeld-Quandt test for Heteroskedasticity',4,TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'Description',header=TRUE)
+ a<-table.element(a,'# significant tests',header=TRUE)
+ a<-table.element(a,'% significant tests',header=TRUE)
+ a<-table.element(a,'OK/NOK',header=TRUE)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'1% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant1)
+ a<-table.element(a,numsignificant1/numgqtests)
+ if (numsignificant1/numgqtests < 0.01) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'5% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant5)
+ a<-table.element(a,numsignificant5/numgqtests)
+ if (numsignificant5/numgqtests < 0.05) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.row.start(a)
+ a<-table.element(a,'10% type I error level',header=TRUE)
+ a<-table.element(a,numsignificant10)
+ a<-table.element(a,numsignificant10/numgqtests)
+ if (numsignificant10/numgqtests < 0.1) dum <- 'OK' else dum <- 'NOK'
+ a<-table.element(a,dum)
+ a<-table.row.end(a)
+ a<-table.end(a)
+ table.save(a,file="/var/www/html/rcomp/tmp/167qho1258725288.tab")
+ }
>
> system("convert tmp/163ay1258725288.ps tmp/163ay1258725288.png")
> system("convert tmp/2nb321258725288.ps tmp/2nb321258725288.png")
> system("convert tmp/3b7l71258725288.ps tmp/3b7l71258725288.png")
> system("convert tmp/43ox61258725288.ps tmp/43ox61258725288.png")
> system("convert tmp/5e5jx1258725288.ps tmp/5e5jx1258725288.png")
> system("convert tmp/6yusr1258725288.ps tmp/6yusr1258725288.png")
> system("convert tmp/75my21258725288.ps tmp/75my21258725288.png")
> system("convert tmp/8kcso1258725288.ps tmp/8kcso1258725288.png")
> system("convert tmp/9ui2d1258725288.ps tmp/9ui2d1258725288.png")
> system("convert tmp/10u5lr1258725288.ps tmp/10u5lr1258725288.png")
>
>
> proc.time()
user system elapsed
2.297 1.574 2.705