## Structural Break
## Using Simulated Data

library(tseries)


set.seed(1000)       
n = 200              

# simulate data, the break date = 100
y = rep(0,n)
e = rnorm(n)
b0 = 1
b1 = 0.2
a0 = 2
a1 = 0.9
br = 100
for (i in 2:br) {
y[i] = b0+b1*y[i-1] + e[i]
}
for (i in (br+1):n) {
y[i] = a0+a1*y[i-1] + e[i]
}

plot(y, type = "l", main = "Series with One Structural Break")

tr = 1:n
d = rep(0,n)
Tb = 100
d[tr>Tb] = 1           
cbind(y,d)

y.1 = c(NA, y[1:n-1])
i.t = d*y.1

# unrestricted model
m.u = lm(y~y.1+d+i.t)
summary(m.u)

# restricted model
m.r = lm(y~y.1)
summary(m.r)

# F test (chow test) with known break date
chowtest = ((0.9899-0.9844)/2)/((1-0.9899)/195)
pvalue = pf(chowtest, 2, 195, lower.tail = FALSE)
c(chowtest, pvalue) 

# grid search the break date
min.br = (0.2*n)
max.br = (0.8*n)

rssv = rep(NA, n)
for (tb in min.br:max.br) {
d = rep(0,n)
d[tr>tb] = 1
y.1 = c(NA, y[1:n-1])
i.t = d*y.1
m.u = lm(y~y.1+d+i.t)
rssv[tb] = sum(resid(m.u)^2)
}

plot(rssv, type = "l", main = "Grid Search: RSS Plot")
min.rss = min(rssv,na.rm = TRUE)
est.breakdate = which(rssv==min.rss) 
est.breakdate