R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i486-pc-linux-gnu (32-bit)
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(10539.51
+ ,2981.85
+ ,11394.84
+ ,10407
+ ,44.23
+ ,10723.78
+ ,3080.58
+ ,11545.71
+ ,10463
+ ,45.85
+ ,10682.06
+ ,3106.22
+ ,11809.38
+ ,10556
+ ,53.38
+ ,10283.19
+ ,3119.31
+ ,11395.64
+ ,10646
+ ,53.26
+ ,10377.18
+ ,3061.26
+ ,11082.38
+ ,10702
+ ,51.8
+ ,10486.64
+ ,3097.31
+ ,11402.75
+ ,11353
+ ,55.3
+ ,10545.38
+ ,3161.69
+ ,11716.87
+ ,11346
+ ,57.81
+ ,10554.27
+ ,3257.16
+ ,12204.98
+ ,11451
+ ,63.96
+ ,10532.54
+ ,3277.01
+ ,12986.62
+ ,11964
+ ,63.77
+ ,10324.31
+ ,3295.32
+ ,13392.79
+ ,12574
+ ,59.15
+ ,10695.25
+ ,3363.99
+ ,14368.05
+ ,13031
+ ,56.12
+ ,10827.81
+ ,3494.17
+ ,15650.83
+ ,13812
+ ,57.42
+ ,10872.48
+ ,3667.03
+ ,16102.64
+ ,14544
+ ,63.52
+ ,10971.19
+ ,3813.06
+ ,16187.64
+ ,14931
+ ,61.71
+ ,11145.65
+ ,3917.96
+ ,16311.54
+ ,14886
+ ,63.01
+ ,11234.68
+ ,3895.51
+ ,17232.97
+ ,16005
+ ,68.18
+ ,11333.88
+ ,3801.06
+ ,16397.83
+ ,17064
+ ,72.03
+ ,10997.97
+ ,3570.12
+ ,14990.31
+ ,15168
+ ,69.75
+ ,11036.89
+ ,3701.61
+ ,15147.55
+ ,16050
+ ,74.41
+ ,11257.35
+ ,3862.27
+ ,15786.78
+ ,15839
+ ,74.33
+ ,11533.59
+ ,3970.1
+ ,15934.09
+ ,15137
+ ,64.24
+ ,11963.12
+ ,4138.52
+ ,16519.44
+ ,14954
+ ,60.03
+ ,12185.15
+ ,4199.75
+ ,16101.07
+ ,15648
+ ,59.44
+ ,12377.62
+ ,4290.89
+ ,16775.08
+ ,15305
+ ,62.5
+ ,12512.89
+ ,4443.91
+ ,17286.32
+ ,15579
+ ,55.04
+ ,12631.48
+ ,4502.64
+ ,17741.23
+ ,16348
+ ,58.34
+ ,12268.53
+ ,4356.98
+ ,17128.37
+ ,15928
+ ,61.92
+ ,12754.8
+ ,4591.27
+ ,17460.53
+ ,16171
+ ,67.65
+ ,13407.75
+ ,4696.96
+ ,17611.14
+ ,15937
+ ,67.68
+ ,13480.21
+ ,4621.4
+ ,18001.37
+ ,15713
+ ,70.3
+ ,13673.28
+ ,4562.84
+ ,17974.77
+ ,15594
+ ,75.26
+ ,13239.71
+ ,4202.52
+ ,16460.95
+ ,15683
+ ,71.44
+ ,13557.69
+ ,4296.49
+ ,16235.39
+ ,16438
+ ,76.36
+ ,13901.28
+ ,4435.23
+ ,16903.36
+ ,17032
+ ,81.71
+ ,13200.58
+ ,4105.18
+ ,15543.76
+ ,17696
+ ,92.6
+ ,13406.97
+ ,4116.68
+ ,15532.18
+ ,17745
+ ,90.6
+ ,12538.12
+ ,3844.49
+ ,13731.31
+ ,19394
+ ,92.23
+ ,12419.57
+ ,3720.98
+ ,13547.84
+ ,20148
+ ,94.09
+ ,12193.88
+ ,3674.4
+ ,12602.93
+ ,20108
+ ,102.79
+ ,12656.63
+ ,3857.62
+ ,13357.7
+ ,18584
+ ,109.65
+ ,12812.48
+ ,3801.06
+ ,13995.33
+ ,18441
+ ,124.05
+ ,12056.67
+ ,3504.37
+ ,14084.6
+ ,18391
+ ,132.69
+ ,11322.38
+ ,3032.6
+ ,13168.91
+ ,19178
+ ,135.81
+ ,11530.75
+ ,3047.03
+ ,12989.35
+ ,18079
+ ,116.07
+ ,11114.08
+ ,2962.34
+ ,12123.53
+ ,18483
+ ,101.42
+ ,9181.73
+ ,2197.82
+ ,9117.03
+ ,19644
+ ,75.73
+ ,8614.55
+ ,2014.45
+ ,8531.45
+ ,19195
+ ,55.48
+ ,8595.56
+ ,1862.83
+ ,8460.94
+ ,19650
+ ,43.8
+ ,8396.2
+ ,1905.41
+ ,8331.49
+ ,20830
+ ,45.29
+ ,7690.5
+ ,1810.99
+ ,7694.78
+ ,23595
+ ,44.01
+ ,7235.47
+ ,1670.07
+ ,7764.58
+ ,22937
+ ,47.48
+ ,7992.12
+ ,1864.44
+ ,8767.96
+ ,21814
+ ,51.07
+ ,8398.37
+ ,2052.02
+ ,9304.43
+ ,21928
+ ,57.84
+ ,8593
+ ,2029.6
+ ,9810.31
+ ,21777
+ ,69.04
+ ,8679.75
+ ,2070.83
+ ,9691.12
+ ,21383
+ ,65.61
+ ,9374.63
+ ,2293.41
+ ,10430.35
+ ,21467
+ ,72.87
+ ,9634.97
+ ,2443.27
+ ,10302.87
+ ,22052
+ ,68.41
+ ,9857.34
+ ,2513.17
+ ,10066.24
+ ,22680
+ ,73.25
+ ,10238.83
+ ,2466.92
+ ,9633.83
+ ,24320
+ ,77.43
+ ,10433.44
+ ,2502.66
+ ,10169.02
+ ,24977
+ ,75.28
+ ,10471.24
+ ,2539.91
+ ,10661.62
+ ,25204
+ ,77.33)
+ ,dim=c(5
+ ,61)
+ ,dimnames=list(c('DowJones'
+ ,'Bel20'
+ ,'Nikkei'
+ ,'Gold'
+ ,'Brent')
+ ,1:61))
> y <- array(NA,dim=c(5,61),dimnames=list(c('DowJones','Bel20','Nikkei','Gold','Brent'),1:61))
> 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 = '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
> 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
DowJones Bel20 Nikkei Gold Brent
1 10539.51 2981.85 11394.84 10407 44.23
2 10723.78 3080.58 11545.71 10463 45.85
3 10682.06 3106.22 11809.38 10556 53.38
4 10283.19 3119.31 11395.64 10646 53.26
5 10377.18 3061.26 11082.38 10702 51.80
6 10486.64 3097.31 11402.75 11353 55.30
7 10545.38 3161.69 11716.87 11346 57.81
8 10554.27 3257.16 12204.98 11451 63.96
9 10532.54 3277.01 12986.62 11964 63.77
10 10324.31 3295.32 13392.79 12574 59.15
11 10695.25 3363.99 14368.05 13031 56.12
12 10827.81 3494.17 15650.83 13812 57.42
13 10872.48 3667.03 16102.64 14544 63.52
14 10971.19 3813.06 16187.64 14931 61.71
15 11145.65 3917.96 16311.54 14886 63.01
16 11234.68 3895.51 17232.97 16005 68.18
17 11333.88 3801.06 16397.83 17064 72.03
18 10997.97 3570.12 14990.31 15168 69.75
19 11036.89 3701.61 15147.55 16050 74.41
20 11257.35 3862.27 15786.78 15839 74.33
21 11533.59 3970.10 15934.09 15137 64.24
22 11963.12 4138.52 16519.44 14954 60.03
23 12185.15 4199.75 16101.07 15648 59.44
24 12377.62 4290.89 16775.08 15305 62.50
25 12512.89 4443.91 17286.32 15579 55.04
26 12631.48 4502.64 17741.23 16348 58.34
27 12268.53 4356.98 17128.37 15928 61.92
28 12754.80 4591.27 17460.53 16171 67.65
29 13407.75 4696.96 17611.14 15937 67.68
30 13480.21 4621.40 18001.37 15713 70.30
31 13673.28 4562.84 17974.77 15594 75.26
32 13239.71 4202.52 16460.95 15683 71.44
33 13557.69 4296.49 16235.39 16438 76.36
34 13901.28 4435.23 16903.36 17032 81.71
35 13200.58 4105.18 15543.76 17696 92.60
36 13406.97 4116.68 15532.18 17745 90.60
37 12538.12 3844.49 13731.31 19394 92.23
38 12419.57 3720.98 13547.84 20148 94.09
39 12193.88 3674.40 12602.93 20108 102.79
40 12656.63 3857.62 13357.70 18584 109.65
41 12812.48 3801.06 13995.33 18441 124.05
42 12056.67 3504.37 14084.60 18391 132.69
43 11322.38 3032.60 13168.91 19178 135.81
44 11530.75 3047.03 12989.35 18079 116.07
45 11114.08 2962.34 12123.53 18483 101.42
46 9181.73 2197.82 9117.03 19644 75.73
47 8614.55 2014.45 8531.45 19195 55.48
48 8595.56 1862.83 8460.94 19650 43.80
49 8396.20 1905.41 8331.49 20830 45.29
50 7690.50 1810.99 7694.78 23595 44.01
51 7235.47 1670.07 7764.58 22937 47.48
52 7992.12 1864.44 8767.96 21814 51.07
53 8398.37 2052.02 9304.43 21928 57.84
54 8593.00 2029.60 9810.31 21777 69.04
55 8679.75 2070.83 9691.12 21383 65.61
56 9374.63 2293.41 10430.35 21467 72.87
57 9634.97 2443.27 10302.87 22052 68.41
58 9857.34 2513.17 10066.24 22680 73.25
59 10238.83 2466.92 9633.83 24320 77.43
60 10433.44 2502.66 10169.02 24977 75.28
61 10471.24 2539.91 10661.62 25204 77.33
> k <- length(x[1,])
> df <- as.data.frame(x)
> (mylm <- lm(df))
Call:
lm(formula = df)
Coefficients:
(Intercept) Bel20 Nikkei Gold Brent
4069.84036 2.69563 -0.27465 0.02823 16.23817
> (mysum <- summary(mylm))
Call:
lm(formula = df)
Residuals:
Min 1Q Median 3Q Max
-622.29 -205.80 -88.22 224.35 759.62
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4069.84036 393.52068 10.342 1.33e-14 ***
Bel20 2.69563 0.21759 12.389 < 2e-16 ***
Nikkei -0.27465 0.05729 -4.794 1.25e-05 ***
Gold 0.02823 0.01565 1.804 0.0767 .
Brent 16.23817 2.63031 6.173 7.93e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 363.2 on 56 degrees of freedom
Multiple R-squared: 0.9541, Adjusted R-squared: 0.9508
F-statistic: 290.8 on 4 and 56 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.0474407734 0.0948815468 0.952559227
[2,] 0.0502882496 0.1005764992 0.949711750
[3,] 0.0348755609 0.0697511217 0.965124439
[4,] 0.0207479977 0.0414959954 0.979252002
[5,] 0.0094231143 0.0188462286 0.990576886
[6,] 0.0037415900 0.0074831800 0.996258410
[7,] 0.0013525637 0.0027051275 0.998647436
[8,] 0.0005105140 0.0010210280 0.999489486
[9,] 0.0008818648 0.0017637297 0.999118135
[10,] 0.0015796625 0.0031593249 0.998420338
[11,] 0.0006786867 0.0013573733 0.999321313
[12,] 0.0003592137 0.0007184275 0.999640786
[13,] 0.0003155094 0.0006310187 0.999684491
[14,] 0.0003231202 0.0006462404 0.999676880
[15,] 0.0004841320 0.0009682640 0.999515868
[16,] 0.0002622009 0.0005244018 0.999737799
[17,] 0.0003753154 0.0007506309 0.999624685
[18,] 0.0002002966 0.0004005932 0.999799703
[19,] 0.0001339745 0.0002679491 0.999866025
[20,] 0.0001555477 0.0003110954 0.999844452
[21,] 0.0006685108 0.0013370217 0.999331489
[22,] 0.0097025105 0.0194050209 0.990297490
[23,] 0.0893604798 0.1787209595 0.910639520
[24,] 0.2756313964 0.5512627928 0.724368604
[25,] 0.6209432976 0.7581134049 0.379056702
[26,] 0.7835915275 0.4328169451 0.216408473
[27,] 0.8235941957 0.3528116086 0.176405804
[28,] 0.7800507459 0.4398985081 0.219949254
[29,] 0.7559702817 0.4880594366 0.244029718
[30,] 0.6929112679 0.6141774642 0.307088732
[31,] 0.6206338750 0.7587322500 0.379366125
[32,] 0.6264587645 0.7470824710 0.373541236
[33,] 0.6844717065 0.6310565870 0.315528294
[34,] 0.7374583687 0.5250832625 0.262541631
[35,] 0.8941786446 0.2116427108 0.105821355
[36,] 0.8497445001 0.3005109998 0.150255500
[37,] 0.8242362015 0.3515275970 0.175763798
[38,] 0.8477495118 0.3045009764 0.152250488
[39,] 0.7841462007 0.4317075985 0.215853799
[40,] 0.7024217401 0.5951565197 0.297578260
[41,] 0.9185719067 0.1628561866 0.081428093
[42,] 0.9978785013 0.0042429974 0.002121499
[43,] 0.9938853285 0.0122293430 0.006114672
[44,] 0.9955555788 0.0088888425 0.004444421
[45,] 0.9895878656 0.0208242689 0.010412134
[46,] 0.9627544881 0.0744910239 0.037245512
> postscript(file="/var/wessaorg/rcomp/tmp/1jjmf1322145782.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/wessaorg/rcomp/tmp/2e4yo1322145782.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/wessaorg/rcomp/tmp/39b421322145782.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/wessaorg/rcomp/tmp/4d7ei1322145782.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/wessaorg/rcomp/tmp/5ham21322145782.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 = 61
Frequency = 1
1 2 3 4 5 6 7
549.27623 480.95644 317.63875 -230.74392 -44.18335 -19.12429 -88.21554
8 9 10 11 12 13 14
-305.44637 -177.40439 -265.63818 224.34870 315.14931 -101.77676 -354.89955
15 16 17 18 19 20 21
-449.02093 -161.94636 -129.93309 -139.33821 -512.25229 -542.05027 -332.35800
22 23 24 25 26 27 28
-122.52918 -190.47227 -98.56867 -121.97027 -112.05039 -296.95253 -450.91891
29 30 31 32 33 34 35
-35.38527 311.71423 578.15306 759.61626 661.12926 710.54169 330.53712
36 37 38 39 40 41 42
533.83980 -168.92412 -56.41829 -556.20981 -448.42018 -194.77140 -265.18281
43 44 45 46 47 48 49
-52.13410 419.59172 219.89852 -92.94927 14.83776 562.00912 154.80501
50 51 52 53 54 55 56
-528.52808 -622.28793 -140.59623 -205.80119 10.60142 20.29580 197.95248
57 58 59 60 61
75.21831 -52.15054 221.07156 482.69323 515.67720
> postscript(file="/var/wessaorg/rcomp/tmp/6zfw91322145782.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 = 61
Frequency = 1
lag(myerror, k = 1) myerror
0 549.27623 NA
1 480.95644 549.27623
2 317.63875 480.95644
3 -230.74392 317.63875
4 -44.18335 -230.74392
5 -19.12429 -44.18335
6 -88.21554 -19.12429
7 -305.44637 -88.21554
8 -177.40439 -305.44637
9 -265.63818 -177.40439
10 224.34870 -265.63818
11 315.14931 224.34870
12 -101.77676 315.14931
13 -354.89955 -101.77676
14 -449.02093 -354.89955
15 -161.94636 -449.02093
16 -129.93309 -161.94636
17 -139.33821 -129.93309
18 -512.25229 -139.33821
19 -542.05027 -512.25229
20 -332.35800 -542.05027
21 -122.52918 -332.35800
22 -190.47227 -122.52918
23 -98.56867 -190.47227
24 -121.97027 -98.56867
25 -112.05039 -121.97027
26 -296.95253 -112.05039
27 -450.91891 -296.95253
28 -35.38527 -450.91891
29 311.71423 -35.38527
30 578.15306 311.71423
31 759.61626 578.15306
32 661.12926 759.61626
33 710.54169 661.12926
34 330.53712 710.54169
35 533.83980 330.53712
36 -168.92412 533.83980
37 -56.41829 -168.92412
38 -556.20981 -56.41829
39 -448.42018 -556.20981
40 -194.77140 -448.42018
41 -265.18281 -194.77140
42 -52.13410 -265.18281
43 419.59172 -52.13410
44 219.89852 419.59172
45 -92.94927 219.89852
46 14.83776 -92.94927
47 562.00912 14.83776
48 154.80501 562.00912
49 -528.52808 154.80501
50 -622.28793 -528.52808
51 -140.59623 -622.28793
52 -205.80119 -140.59623
53 10.60142 -205.80119
54 20.29580 10.60142
55 197.95248 20.29580
56 75.21831 197.95248
57 -52.15054 75.21831
58 221.07156 -52.15054
59 482.69323 221.07156
60 515.67720 482.69323
61 NA 515.67720
> dum1 <- dum[2:length(myerror),]
> dum1
lag(myerror, k = 1) myerror
[1,] 480.95644 549.27623
[2,] 317.63875 480.95644
[3,] -230.74392 317.63875
[4,] -44.18335 -230.74392
[5,] -19.12429 -44.18335
[6,] -88.21554 -19.12429
[7,] -305.44637 -88.21554
[8,] -177.40439 -305.44637
[9,] -265.63818 -177.40439
[10,] 224.34870 -265.63818
[11,] 315.14931 224.34870
[12,] -101.77676 315.14931
[13,] -354.89955 -101.77676
[14,] -449.02093 -354.89955
[15,] -161.94636 -449.02093
[16,] -129.93309 -161.94636
[17,] -139.33821 -129.93309
[18,] -512.25229 -139.33821
[19,] -542.05027 -512.25229
[20,] -332.35800 -542.05027
[21,] -122.52918 -332.35800
[22,] -190.47227 -122.52918
[23,] -98.56867 -190.47227
[24,] -121.97027 -98.56867
[25,] -112.05039 -121.97027
[26,] -296.95253 -112.05039
[27,] -450.91891 -296.95253
[28,] -35.38527 -450.91891
[29,] 311.71423 -35.38527
[30,] 578.15306 311.71423
[31,] 759.61626 578.15306
[32,] 661.12926 759.61626
[33,] 710.54169 661.12926
[34,] 330.53712 710.54169
[35,] 533.83980 330.53712
[36,] -168.92412 533.83980
[37,] -56.41829 -168.92412
[38,] -556.20981 -56.41829
[39,] -448.42018 -556.20981
[40,] -194.77140 -448.42018
[41,] -265.18281 -194.77140
[42,] -52.13410 -265.18281
[43,] 419.59172 -52.13410
[44,] 219.89852 419.59172
[45,] -92.94927 219.89852
[46,] 14.83776 -92.94927
[47,] 562.00912 14.83776
[48,] 154.80501 562.00912
[49,] -528.52808 154.80501
[50,] -622.28793 -528.52808
[51,] -140.59623 -622.28793
[52,] -205.80119 -140.59623
[53,] 10.60142 -205.80119
[54,] 20.29580 10.60142
[55,] 197.95248 20.29580
[56,] 75.21831 197.95248
[57,] -52.15054 75.21831
[58,] 221.07156 -52.15054
[59,] 482.69323 221.07156
[60,] 515.67720 482.69323
> z <- as.data.frame(dum1)
> z
lag(myerror, k = 1) myerror
1 480.95644 549.27623
2 317.63875 480.95644
3 -230.74392 317.63875
4 -44.18335 -230.74392
5 -19.12429 -44.18335
6 -88.21554 -19.12429
7 -305.44637 -88.21554
8 -177.40439 -305.44637
9 -265.63818 -177.40439
10 224.34870 -265.63818
11 315.14931 224.34870
12 -101.77676 315.14931
13 -354.89955 -101.77676
14 -449.02093 -354.89955
15 -161.94636 -449.02093
16 -129.93309 -161.94636
17 -139.33821 -129.93309
18 -512.25229 -139.33821
19 -542.05027 -512.25229
20 -332.35800 -542.05027
21 -122.52918 -332.35800
22 -190.47227 -122.52918
23 -98.56867 -190.47227
24 -121.97027 -98.56867
25 -112.05039 -121.97027
26 -296.95253 -112.05039
27 -450.91891 -296.95253
28 -35.38527 -450.91891
29 311.71423 -35.38527
30 578.15306 311.71423
31 759.61626 578.15306
32 661.12926 759.61626
33 710.54169 661.12926
34 330.53712 710.54169
35 533.83980 330.53712
36 -168.92412 533.83980
37 -56.41829 -168.92412
38 -556.20981 -56.41829
39 -448.42018 -556.20981
40 -194.77140 -448.42018
41 -265.18281 -194.77140
42 -52.13410 -265.18281
43 419.59172 -52.13410
44 219.89852 419.59172
45 -92.94927 219.89852
46 14.83776 -92.94927
47 562.00912 14.83776
48 154.80501 562.00912
49 -528.52808 154.80501
50 -622.28793 -528.52808
51 -140.59623 -622.28793
52 -205.80119 -140.59623
53 10.60142 -205.80119
54 20.29580 10.60142
55 197.95248 20.29580
56 75.21831 197.95248
57 -52.15054 75.21831
58 221.07156 -52.15054
59 482.69323 221.07156
60 515.67720 482.69323
> 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/wessaorg/rcomp/tmp/757tb1322145782.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/wessaorg/rcomp/tmp/8ftey1322145782.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/wessaorg/rcomp/tmp/9d0dr1322145782.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/wessaorg/rcomp/tmp/105cge1322145782.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/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab
> load(file="/var/wessaorg/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/wessaorg/rcomp/tmp/11nkg11322145782.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/wessaorg/rcomp/tmp/12svka1322145782.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/wessaorg/rcomp/tmp/1356c41322145783.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/wessaorg/rcomp/tmp/14m8dr1322145783.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/wessaorg/rcomp/tmp/15g72w1322145783.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/wessaorg/rcomp/tmp/161hvf1322145783.tab")
+ }
>
> try(system("convert tmp/1jjmf1322145782.ps tmp/1jjmf1322145782.png",intern=TRUE))
character(0)
> try(system("convert tmp/2e4yo1322145782.ps tmp/2e4yo1322145782.png",intern=TRUE))
character(0)
> try(system("convert tmp/39b421322145782.ps tmp/39b421322145782.png",intern=TRUE))
character(0)
> try(system("convert tmp/4d7ei1322145782.ps tmp/4d7ei1322145782.png",intern=TRUE))
character(0)
> try(system("convert tmp/5ham21322145782.ps tmp/5ham21322145782.png",intern=TRUE))
character(0)
> try(system("convert tmp/6zfw91322145782.ps tmp/6zfw91322145782.png",intern=TRUE))
character(0)
> try(system("convert tmp/757tb1322145782.ps tmp/757tb1322145782.png",intern=TRUE))
character(0)
> try(system("convert tmp/8ftey1322145782.ps tmp/8ftey1322145782.png",intern=TRUE))
character(0)
> try(system("convert tmp/9d0dr1322145782.ps tmp/9d0dr1322145782.png",intern=TRUE))
character(0)
> try(system("convert tmp/105cge1322145782.ps tmp/105cge1322145782.png",intern=TRUE))
character(0)
>
>
> proc.time()
user system elapsed
3.494 0.545 4.110