# Code to simulate HMM-GARCH(1,1) Process
objects();
remove(list=objects());
options(scipen=20);

# Begin with simulating the hidden markov process.
n = 500 #Length of the data set

# Initial distribution in the chain
v1 = 0.5
v2 = 0.5

# Transition probabilities
p12 = 0.01
p21 = 0.01

# Parameter values
w1 <- 0.07
w2 <- 7.75
a1 <- 0.04
b1 <- 0.03
a2 <- 0.21
b2 <- 0.77

startState = sample(1,c(1,2))
transProbs =c(0,runif(n-1)) # transition probabilities
hmp <- rep(0,n)
hmp[1] = startState
for(i in 2:n){
	prevState = hmp[i-1]
	if (prevState == 1 & transProbs[i] <= p12){
		hmp[i]=2} #switch to state 2
	else if (prevState == 2 & transProbs[i] <= p21){
		hmp[i]=1}	#switch to state 1
	else{
		hmp[i] = hmp[i-1]} #stay in the same state
}

# Generate the HMM training data
# Generate the squared innovations
eps = rnorm(n)
h <- rep(0,n)

h[1] = 0.5*(w1/(1-a1-b1) + w2/(1-a2-b2)) 
rt <- rep(0,n) # returns
rt[1] = eps[1]*sqrt(h[1])

for (i in 2:n){
	prevState = hmp[i-1]
	if (prevState == 1){
		h[i] = w1 + b1*h[i-1] + a1*h[i-1]*(eps[i-1])^2
	}
	if (prevState == 2){
		h[i] = w2 + b2*h[i-1] + a2*h[i-1]*(eps[i-1])^2

	}
  	rt[i] = eps[i]*sqrt(h[i])
}


write.table(rt,"returns.txt",sep=" ",row.names=FALSE,col.names = FALSE)
write.table(hmp,'hmp.txt',sep=' ',row.names=FALSE,col.names = FALSE)
jpeg('blah.jpeg',width=1600,height=800)
par(mfrow=c(3,1))
plot(hmp,type='l')
plot(rt,type='l')
plot(h, type = 'l')
dev.off()