[R] Break Points
Ayanendu Sanyal
ayanendu1 at gmail.com
Mon Jan 2 13:21:30 CET 2012
Respected Sir
I fitted that model... lnrpe= a+bt and conducted sup F test since
autocorelation is present in the data and AIC as you mentioned might
not .... Is it OK... Since I am not well versed with time series
econometrics can you please tell me if the
work is now correct or not
--
Please have a look at our new mission and contribute into it (cut and
paste the link below in the address bar of your internet browser)
http://thesocialscienceinformer.blogspot.com/
Thanking you
Ayanendu Sanyal
PhD Scholar
Institute for Social and Economic Change (ISEC)
P.O- Nagarbhavi
Bangalore-72
State- Karnataka
Country- India
PIN- 560072
www.isec.ac.in/phd.html
http://ayanendusanyal.blogspot.com/
-------------- next part --------------
lnrpe
1.6113515
1.619601724
1.599889264
1.645954835
1.723777317
1.830606002
2.034751407
2.112377045
2.095050993
2.046822835
2.276628064
2.543864584
2.619425807
2.717786454
2.874537082
2.923223972
3.136825311
3.206377996
3.352655132
3.49032806
3.508602739
3.621768106
3.803617305
4.141727497
4.27471221
4.34523451
4.242555261
4.262046942
4.378894917
4.419018243
4.391496862
4.489015146
4.588180885
4.846861554
5.069645376
5.257766481
5.292695491
5.307844982
5.277006289
5.323613873
5.377629069
5.429256311
5.443411354
5.47242567
5.727392687
6.147054773
-------------- next part --------------
R version 2.14.1 (2011-12-22)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i386-pc-mingw32/i386 (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.
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.
> source("C:\\Documents and Settings\\ayaendu\\My Documents\\struc.R")
Loading required package: MASS
Loading required package: nnet
Loading required package: survival
Loading required package: splines
Loading required package: zoo
Attaching package: ?zoo?
The following object(s) are masked from ?package:base?:
as.Date, as.Date.numeric
Loading required package: sandwich
Warning message:
Overlapping confidence intervals
> ayan <- read.table("ayan.txt", header=T) ## Fetching files saved in My documents
> ayan = ts(ayan, start=1, frequency=1) ## Setting the data into time series data
> trend = time (ayan) ## defining the explanatory variable
> reg = lm(ayan~trend) ## regression lnrpe= a+bt
> summary (reg) ## regression results summary
Call:
lm(formula = ayan ~ trend)
Residuals:
Min 1Q Median 3Q Max
-0.35086 -0.11196 -0.01008 0.09193 0.38507
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.345373 0.049774 27.03 <2e-16 ***
trend 0.101771 0.001844 55.19 <2e-16 ***
---
Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
Residual standard error: 0.166 on 44 degrees of freedom
Multiple R-squared: 0.9858, Adjusted R-squared: 0.9854
F-statistic: 3046 on 1 and 44 DF, p-value: < 2.2e-16
> library (car) ## loading package to do dwtest
> library (strucchange) ## loading package to do structural break
> library (lmtest) ## loading package to do dwtest
> dwtest (reg) ## for testing serial auto correrelation
Durbin-Watson test
data: reg
DW = 0.4157, p-value = 1.405e-12
alternative hypothesis: true autocorrelation is greater than 0
> bp.ayan <- breakpoints(ayan ~ trend) ## searching for break points
> summary (bp.ayan) ## checking the breakpoints and minimum AIC level
Optimal (m+1)-segment partition:
Call:
breakpoints.formula(formula = ayan ~ trend)
Breakpoints at observation number:
m = 1 23
m = 2 23 33
m = 3 11 23 33
m = 4 11 23 33 40
m = 5 11 20 26 34 40
m = 6 8 14 20 26 34 40
Corresponding to breakdates:
m = 1 23
m = 2 23 33
m = 3 11 23 33
m = 4 11 23 33 40
m = 5 11 20 26 34 40
m = 6 8 14 20 26 34 40
Fit:
m 0 1 2 3 4 5
RSS 1.2131579 0.7466773 0.5243361 0.3570253 0.2712234 0.2589809
BIC -25.2008012 -36.0409283 -40.8160181 -47.0090984 -48.1669059 -38.8056613
m 6
RSS 0.2698858
BIC -25.4224766
> ci.ayan <- confint(bp.ayan)## confidence interval for the test
> breakdates (ci.ayan)
2.5 % breakpoints 97.5 %
1 10 11 15
2 22 23 24
3 32 33 34
4 25 40 41
Warning message:
Overlapping confidence intervals
> ci.ayan ## displays the confidence interval
Confidence intervals for breakpoints
of optimal 5-segment partition:
Call:
confint.breakpointsfull(object = bp.ayan)
Breakpoints at observation number:
2.5 % breakpoints 97.5 %
1 10 11 15
2 22 23 24
3 32 33 34
4 25 40 41
Corresponding to breakdates:
2.5 % breakpoints 97.5 %
1 10 11 15
2 22 23 24
3 32 33 34
4 25 40 41
Warning message:
Overlapping confidence intervals
> fs.ayan <- Fstats(ayan~trend, data = ayan, from = 0.1) ##sup F test h=0.1 trimming parameter
> plot (fs.ayan)
> bp3 <- breakpoints(ayan ~ trend,breaks=3)## searching for break points with breakpoints assumed as three since AIC min can overestimate in presence of serial auto corelation
> summary (bp3)## checking the break points
Optimal (m+1)-segment partition:
Call:
breakpoints.formula(formula = ayan ~ trend, breaks = 3)
Breakpoints at observation number:
m = 1 23
m = 2 23 33
m = 3 11 23 33
Corresponding to breakdates:
m = 1 23
m = 2 23 33
m = 3 11 23 33
Fit:
m 0 1 2 3
RSS 1.2131579 0.7466773 0.5243361 0.3570253
BIC -25.2008012 -36.0409283 -40.8160181 -47.0090984
> ci.ayan1 <- confint(bp3)## confidence interval for the test when break points are 3
> breakdates (ci.ayan1) ## displays the confidence interval
2.5 % breakpoints 97.5 %
1 10 11 15
2 22 23 24
3 32 33 34
> ci.ayan1 ## displays the confidence interval for three breaks
Confidence intervals for breakpoints
of optimal 4-segment partition:
Call:
confint.breakpointsfull(object = bp3)
Breakpoints at observation number:
2.5 % breakpoints 97.5 %
1 10 11 15
2 22 23 24
3 32 33 34
Corresponding to breakdates:
2.5 % breakpoints 97.5 %
1 10 11 15
2 22 23 24
3 32 33 34
> plot (ayan) ## plots the series
> lines (ci.ayan1) ## shows the breaks and the CIs for three breaks
>
>
More information about the R-help
mailing list