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(87032
+ ,537000
+ ,88219
+ ,90390
+ ,90269
+ ,90398
+ ,87175
+ ,543000
+ ,87032
+ ,88219
+ ,90390
+ ,90269
+ ,92603
+ ,594000
+ ,87175
+ ,87032
+ ,88219
+ ,90390
+ ,93571
+ ,611000
+ ,92603
+ ,87175
+ ,87032
+ ,88219
+ ,94118
+ ,613000
+ ,93571
+ ,92603
+ ,87175
+ ,87032
+ ,92159
+ ,611000
+ ,94118
+ ,93571
+ ,92603
+ ,87175
+ ,89528
+ ,594000
+ ,92159
+ ,94118
+ ,93571
+ ,92603
+ ,89955
+ ,595000
+ ,89528
+ ,92159
+ ,94118
+ ,93571
+ ,89587
+ ,591000
+ ,89955
+ ,89528
+ ,92159
+ ,94118
+ ,89488
+ ,589000
+ ,89587
+ ,89955
+ ,89528
+ ,92159
+ ,88521
+ ,584000
+ ,89488
+ ,89587
+ ,89955
+ ,89528
+ ,86587
+ ,573000
+ ,88521
+ ,89488
+ ,89587
+ ,89955
+ ,85159
+ ,567000
+ ,86587
+ ,88521
+ ,89488
+ ,89587
+ ,84915
+ ,569000
+ ,85159
+ ,86587
+ ,88521
+ ,89488
+ ,91378
+ ,621000
+ ,84915
+ ,85159
+ ,86587
+ ,88521
+ ,92729
+ ,629000
+ ,91378
+ ,84915
+ ,85159
+ ,86587
+ ,92194
+ ,628000
+ ,92729
+ ,91378
+ ,84915
+ ,85159
+ ,89664
+ ,612000
+ ,92194
+ ,92729
+ ,91378
+ ,84915
+ ,86285
+ ,595000
+ ,89664
+ ,92194
+ ,92729
+ ,91378
+ ,86858
+ ,597000
+ ,86285
+ ,89664
+ ,92194
+ ,92729
+ ,87184
+ ,593000
+ ,86858
+ ,86285
+ ,89664
+ ,92194
+ ,86629
+ ,590000
+ ,87184
+ ,86858
+ ,86285
+ ,89664
+ ,85220
+ ,580000
+ ,86629
+ ,87184
+ ,86858
+ ,86285
+ ,84816
+ ,574000
+ ,85220
+ ,86629
+ ,87184
+ ,86858
+ ,84831
+ ,573000
+ ,84816
+ ,85220
+ ,86629
+ ,87184
+ ,84957
+ ,573000
+ ,84831
+ ,84816
+ ,85220
+ ,86629
+ ,90951
+ ,620000
+ ,84957
+ ,84831
+ ,84816
+ ,85220
+ ,92134
+ ,626000
+ ,90951
+ ,84957
+ ,84831
+ ,84816
+ ,91790
+ ,620000
+ ,92134
+ ,90951
+ ,84957
+ ,84831
+ ,86625
+ ,588000
+ ,91790
+ ,92134
+ ,90951
+ ,84957
+ ,83324
+ ,566000
+ ,86625
+ ,91790
+ ,92134
+ ,90951
+ ,82719
+ ,557000
+ ,83324
+ ,86625
+ ,91790
+ ,92134
+ ,83614
+ ,561000
+ ,82719
+ ,83324
+ ,86625
+ ,91790
+ ,81640
+ ,549000
+ ,83614
+ ,82719
+ ,83324
+ ,86625
+ ,78665
+ ,532000
+ ,81640
+ ,83614
+ ,82719
+ ,83324
+ ,77828
+ ,526000
+ ,78665
+ ,81640
+ ,83614
+ ,82719
+ ,75728
+ ,511000
+ ,77828
+ ,78665
+ ,81640
+ ,83614
+ ,72187
+ ,499000
+ ,75728
+ ,77828
+ ,78665
+ ,81640
+ ,79357
+ ,555000
+ ,72187
+ ,75728
+ ,77828
+ ,78665
+ ,81329
+ ,565000
+ ,79357
+ ,72187
+ ,75728
+ ,77828
+ ,77304
+ ,542000
+ ,81329
+ ,79357
+ ,72187
+ ,75728
+ ,75576
+ ,527000
+ ,77304
+ ,81329
+ ,79357
+ ,72187
+ ,72932
+ ,510000
+ ,75576
+ ,77304
+ ,81329
+ ,79357
+ ,74291
+ ,514000
+ ,72932
+ ,75576
+ ,77304
+ ,81329
+ ,74988
+ ,517000
+ ,74291
+ ,72932
+ ,75576
+ ,77304
+ ,73302
+ ,508000
+ ,74988
+ ,74291
+ ,72932
+ ,75576
+ ,70483
+ ,493000
+ ,73302
+ ,74988
+ ,74291
+ ,72932
+ ,69848
+ ,490000
+ ,70483
+ ,73302
+ ,74988
+ ,74291
+ ,66466
+ ,469000
+ ,69848
+ ,70483
+ ,73302
+ ,74988
+ ,67610
+ ,478000
+ ,66466
+ ,69848
+ ,70483
+ ,73302
+ ,75091
+ ,528000
+ ,67610
+ ,66466
+ ,69848
+ ,70483
+ ,76207
+ ,534000
+ ,75091
+ ,67610
+ ,66466
+ ,69848
+ ,73454
+ ,518000
+ ,76207
+ ,75091
+ ,67610
+ ,66466
+ ,72008
+ ,506000
+ ,73454
+ ,76207
+ ,75091
+ ,67610
+ ,71362
+ ,502000
+ ,72008
+ ,73454
+ ,76207
+ ,75091
+ ,74250
+ ,516000
+ ,71362
+ ,72008
+ ,73454
+ ,76207)
+ ,dim=c(6
+ ,56)
+ ,dimnames=list(c('Y'
+ ,'X'
+ ,'Y1'
+ ,'Y2'
+ ,'Y3'
+ ,'Y4')
+ ,1:56))
> y <- array(NA,dim=c(6,56),dimnames=list(c('Y','X','Y1','Y2','Y3','Y4'),1:56))
> 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 = 'Do not include Seasonal 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 t
1 87032 537000 88219 90390 90269 90398 1
2 87175 543000 87032 88219 90390 90269 2
3 92603 594000 87175 87032 88219 90390 3
4 93571 611000 92603 87175 87032 88219 4
5 94118 613000 93571 92603 87175 87032 5
6 92159 611000 94118 93571 92603 87175 6
7 89528 594000 92159 94118 93571 92603 7
8 89955 595000 89528 92159 94118 93571 8
9 89587 591000 89955 89528 92159 94118 9
10 89488 589000 89587 89955 89528 92159 10
11 88521 584000 89488 89587 89955 89528 11
12 86587 573000 88521 89488 89587 89955 12
13 85159 567000 86587 88521 89488 89587 13
14 84915 569000 85159 86587 88521 89488 14
15 91378 621000 84915 85159 86587 88521 15
16 92729 629000 91378 84915 85159 86587 16
17 92194 628000 92729 91378 84915 85159 17
18 89664 612000 92194 92729 91378 84915 18
19 86285 595000 89664 92194 92729 91378 19
20 86858 597000 86285 89664 92194 92729 20
21 87184 593000 86858 86285 89664 92194 21
22 86629 590000 87184 86858 86285 89664 22
23 85220 580000 86629 87184 86858 86285 23
24 84816 574000 85220 86629 87184 86858 24
25 84831 573000 84816 85220 86629 87184 25
26 84957 573000 84831 84816 85220 86629 26
27 90951 620000 84957 84831 84816 85220 27
28 92134 626000 90951 84957 84831 84816 28
29 91790 620000 92134 90951 84957 84831 29
30 86625 588000 91790 92134 90951 84957 30
31 83324 566000 86625 91790 92134 90951 31
32 82719 557000 83324 86625 91790 92134 32
33 83614 561000 82719 83324 86625 91790 33
34 81640 549000 83614 82719 83324 86625 34
35 78665 532000 81640 83614 82719 83324 35
36 77828 526000 78665 81640 83614 82719 36
37 75728 511000 77828 78665 81640 83614 37
38 72187 499000 75728 77828 78665 81640 38
39 79357 555000 72187 75728 77828 78665 39
40 81329 565000 79357 72187 75728 77828 40
41 77304 542000 81329 79357 72187 75728 41
42 75576 527000 77304 81329 79357 72187 42
43 72932 510000 75576 77304 81329 79357 43
44 74291 514000 72932 75576 77304 81329 44
45 74988 517000 74291 72932 75576 77304 45
46 73302 508000 74988 74291 72932 75576 46
47 70483 493000 73302 74988 74291 72932 47
48 69848 490000 70483 73302 74988 74291 48
49 66466 469000 69848 70483 73302 74988 49
50 67610 478000 66466 69848 70483 73302 50
51 75091 528000 67610 66466 69848 70483 51
52 76207 534000 75091 67610 66466 69848 52
53 73454 518000 76207 75091 67610 66466 53
54 72008 506000 73454 76207 75091 67610 54
55 71362 502000 72008 73454 76207 75091 55
56 74250 516000 71362 72008 73454 76207 56
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) X Y1 Y2 Y3 Y4
1.666e+04 1.097e-01 1.987e-01 -1.203e-01 -1.636e-02 4.562e-02
t
-1.544e+02
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-2089.4720 -756.9685 0.4767 600.9534 2316.3675
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.666e+04 4.582e+03 3.636 0.000663 ***
X 1.097e-01 7.528e-03 14.572 < 2e-16 ***
Y1 1.987e-01 8.506e-02 2.336 0.023630 *
Y2 -1.203e-01 9.357e-02 -1.286 0.204550
Y3 -1.636e-02 9.011e-02 -0.182 0.856695
Y4 4.562e-02 6.063e-02 0.752 0.455352
t -1.544e+02 2.100e+01 -7.349 1.91e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1035 on 49 degrees of freedom
Multiple R-squared: 0.9841, Adjusted R-squared: 0.9821
F-statistic: 504.9 on 6 and 49 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.10688205 0.213764093 0.893117954
[2,] 0.05668408 0.113368150 0.943315925
[3,] 0.07816009 0.156320182 0.921839909
[4,] 0.12936512 0.258730241 0.870634879
[5,] 0.13395354 0.267907089 0.866046456
[6,] 0.10569410 0.211388200 0.894305900
[7,] 0.17056880 0.341137603 0.829431198
[8,] 0.12890660 0.257813210 0.871093395
[9,] 0.23387456 0.467749111 0.766125444
[10,] 0.29606125 0.592122504 0.703938748
[11,] 0.44908752 0.898175035 0.550912482
[12,] 0.61505912 0.769881763 0.384940882
[13,] 0.60931057 0.781378851 0.390689426
[14,] 0.56371351 0.872572983 0.436286492
[15,] 0.57264675 0.854706504 0.427353252
[16,] 0.58486411 0.830271773 0.415135886
[17,] 0.58849161 0.823016785 0.411508393
[18,] 0.82468134 0.350637327 0.175318664
[19,] 0.92150586 0.156988284 0.078494142
[20,] 0.97208599 0.055828019 0.027914010
[21,] 0.96976343 0.060473146 0.030236573
[22,] 0.98076191 0.038476183 0.019238092
[23,] 0.98371430 0.032571392 0.016285696
[24,] 0.97512668 0.049746634 0.024873317
[25,] 0.96430832 0.071383363 0.035691681
[26,] 0.96092212 0.078155752 0.039077876
[27,] 0.97427011 0.051459777 0.025729889
[28,] 0.99726166 0.005476677 0.002738339
[29,] 0.99766534 0.004669315 0.002334658
[30,] 0.99595938 0.008081232 0.004040616
[31,] 0.99185763 0.016284738 0.008142369
[32,] 0.99735881 0.005282379 0.002641189
[33,] 0.99243719 0.015125623 0.007562812
[34,] 0.99767037 0.004659250 0.002329625
[35,] 0.99779886 0.004402284 0.002201142
[36,] 0.99079485 0.018410307 0.009205154
[37,] 0.96450762 0.070984754 0.035492377
> postscript(file="/var/www/html/rcomp/tmp/1g9781258572573.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/2biv41258572573.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/39c6r1258572573.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/43qnj1258572573.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/5zh8k1258572573.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 = 56
Frequency = 1
1 2 3 4 5 6
2316.36754 1938.02700 1713.36502 -10.86891 988.30269 -506.89183
7 8 9 10 11 12
-895.37821 -171.87463 -405.10646 40.47747 -121.24855 -539.46007
13 14 15 16 17 18
-871.80201 -1141.10620 -339.05546 -959.94320 -660.57901 -895.30538
19 20 21 22 23 24
-2089.47200 -1284.92031 -903.13212 -910.35961 -754.96429 -154.02983
25 26 27 28 29 30
11.82231 242.85975 1269.70631 791.72567 1747.76691 550.54177
31 32 33 34 35 36
548.06467 1059.59365 1324.39805 752.18840 437.06509 808.48014
37 38 39 40 41 42
243.59756 -1468.69125 285.35181 -532.12893 -1370.94411 16.76644
43 44 45 46 47 48
-1043.74221 192.43499 281.91881 -201.81710 -659.23382 -504.11108
49 50 51 52 53 54
-1700.40951 -762.98098 871.33872 108.34896 116.25509 892.47934
55 56
472.67199 1837.64090
> postscript(file="/var/www/html/rcomp/tmp/6w5h71258572573.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 = 56
Frequency = 1
lag(myerror, k = 1) myerror
0 2316.36754 NA
1 1938.02700 2316.36754
2 1713.36502 1938.02700
3 -10.86891 1713.36502
4 988.30269 -10.86891
5 -506.89183 988.30269
6 -895.37821 -506.89183
7 -171.87463 -895.37821
8 -405.10646 -171.87463
9 40.47747 -405.10646
10 -121.24855 40.47747
11 -539.46007 -121.24855
12 -871.80201 -539.46007
13 -1141.10620 -871.80201
14 -339.05546 -1141.10620
15 -959.94320 -339.05546
16 -660.57901 -959.94320
17 -895.30538 -660.57901
18 -2089.47200 -895.30538
19 -1284.92031 -2089.47200
20 -903.13212 -1284.92031
21 -910.35961 -903.13212
22 -754.96429 -910.35961
23 -154.02983 -754.96429
24 11.82231 -154.02983
25 242.85975 11.82231
26 1269.70631 242.85975
27 791.72567 1269.70631
28 1747.76691 791.72567
29 550.54177 1747.76691
30 548.06467 550.54177
31 1059.59365 548.06467
32 1324.39805 1059.59365
33 752.18840 1324.39805
34 437.06509 752.18840
35 808.48014 437.06509
36 243.59756 808.48014
37 -1468.69125 243.59756
38 285.35181 -1468.69125
39 -532.12893 285.35181
40 -1370.94411 -532.12893
41 16.76644 -1370.94411
42 -1043.74221 16.76644
43 192.43499 -1043.74221
44 281.91881 192.43499
45 -201.81710 281.91881
46 -659.23382 -201.81710
47 -504.11108 -659.23382
48 -1700.40951 -504.11108
49 -762.98098 -1700.40951
50 871.33872 -762.98098
51 108.34896 871.33872
52 116.25509 108.34896
53 892.47934 116.25509
54 472.67199 892.47934
55 1837.64090 472.67199
56 NA 1837.64090
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 1938.02700 2316.36754
[2,] 1713.36502 1938.02700
[3,] -10.86891 1713.36502
[4,] 988.30269 -10.86891
[5,] -506.89183 988.30269
[6,] -895.37821 -506.89183
[7,] -171.87463 -895.37821
[8,] -405.10646 -171.87463
[9,] 40.47747 -405.10646
[10,] -121.24855 40.47747
[11,] -539.46007 -121.24855
[12,] -871.80201 -539.46007
[13,] -1141.10620 -871.80201
[14,] -339.05546 -1141.10620
[15,] -959.94320 -339.05546
[16,] -660.57901 -959.94320
[17,] -895.30538 -660.57901
[18,] -2089.47200 -895.30538
[19,] -1284.92031 -2089.47200
[20,] -903.13212 -1284.92031
[21,] -910.35961 -903.13212
[22,] -754.96429 -910.35961
[23,] -154.02983 -754.96429
[24,] 11.82231 -154.02983
[25,] 242.85975 11.82231
[26,] 1269.70631 242.85975
[27,] 791.72567 1269.70631
[28,] 1747.76691 791.72567
[29,] 550.54177 1747.76691
[30,] 548.06467 550.54177
[31,] 1059.59365 548.06467
[32,] 1324.39805 1059.59365
[33,] 752.18840 1324.39805
[34,] 437.06509 752.18840
[35,] 808.48014 437.06509
[36,] 243.59756 808.48014
[37,] -1468.69125 243.59756
[38,] 285.35181 -1468.69125
[39,] -532.12893 285.35181
[40,] -1370.94411 -532.12893
[41,] 16.76644 -1370.94411
[42,] -1043.74221 16.76644
[43,] 192.43499 -1043.74221
[44,] 281.91881 192.43499
[45,] -201.81710 281.91881
[46,] -659.23382 -201.81710
[47,] -504.11108 -659.23382
[48,] -1700.40951 -504.11108
[49,] -762.98098 -1700.40951
[50,] 871.33872 -762.98098
[51,] 108.34896 871.33872
[52,] 116.25509 108.34896
[53,] 892.47934 116.25509
[54,] 472.67199 892.47934
[55,] 1837.64090 472.67199
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 1938.02700 2316.36754
2 1713.36502 1938.02700
3 -10.86891 1713.36502
4 988.30269 -10.86891
5 -506.89183 988.30269
6 -895.37821 -506.89183
7 -171.87463 -895.37821
8 -405.10646 -171.87463
9 40.47747 -405.10646
10 -121.24855 40.47747
11 -539.46007 -121.24855
12 -871.80201 -539.46007
13 -1141.10620 -871.80201
14 -339.05546 -1141.10620
15 -959.94320 -339.05546
16 -660.57901 -959.94320
17 -895.30538 -660.57901
18 -2089.47200 -895.30538
19 -1284.92031 -2089.47200
20 -903.13212 -1284.92031
21 -910.35961 -903.13212
22 -754.96429 -910.35961
23 -154.02983 -754.96429
24 11.82231 -154.02983
25 242.85975 11.82231
26 1269.70631 242.85975
27 791.72567 1269.70631
28 1747.76691 791.72567
29 550.54177 1747.76691
30 548.06467 550.54177
31 1059.59365 548.06467
32 1324.39805 1059.59365
33 752.18840 1324.39805
34 437.06509 752.18840
35 808.48014 437.06509
36 243.59756 808.48014
37 -1468.69125 243.59756
38 285.35181 -1468.69125
39 -532.12893 285.35181
40 -1370.94411 -532.12893
41 16.76644 -1370.94411
42 -1043.74221 16.76644
43 192.43499 -1043.74221
44 281.91881 192.43499
45 -201.81710 281.91881
46 -659.23382 -201.81710
47 -504.11108 -659.23382
48 -1700.40951 -504.11108
49 -762.98098 -1700.40951
50 871.33872 -762.98098
51 108.34896 871.33872
52 116.25509 108.34896
53 892.47934 116.25509
54 472.67199 892.47934
55 1837.64090 472.67199
> 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/7iyir1258572573.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/8gk1d1258572573.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/9wa011258572573.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/10s7gw1258572573.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/11sg4b1258572573.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/12pzlg1258572573.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/133xzk1258572573.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/14uj9p1258572573.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/15e2fi1258572573.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/16yelr1258572573.tab")
+ }
>
> system("convert tmp/1g9781258572573.ps tmp/1g9781258572573.png")
> system("convert tmp/2biv41258572573.ps tmp/2biv41258572573.png")
> system("convert tmp/39c6r1258572573.ps tmp/39c6r1258572573.png")
> system("convert tmp/43qnj1258572573.ps tmp/43qnj1258572573.png")
> system("convert tmp/5zh8k1258572573.ps tmp/5zh8k1258572573.png")
> system("convert tmp/6w5h71258572573.ps tmp/6w5h71258572573.png")
> system("convert tmp/7iyir1258572573.ps tmp/7iyir1258572573.png")
> system("convert tmp/8gk1d1258572573.ps tmp/8gk1d1258572573.png")
> system("convert tmp/9wa011258572573.ps tmp/9wa011258572573.png")
> system("convert tmp/10s7gw1258572573.ps tmp/10s7gw1258572573.png")
>
>
> proc.time()
user system elapsed
2.462 1.564 3.229