Hamilton's "basic filter" algorithm implemented in R

GtownMath611hw11prob1bBasicFilter.R

index2slag<-function(index){
  # Takes an index between 1 and 32 and converts it into a binary vector. This is used to iterate all of the possible values of (s_t,s_t1,s_t2,s_t3,s_t4)
  num<-paste(rev(intToBits(index-1)),"")
  return(as.logical(as.numeric(num[28:32])))
}

ref.s<-t(sapply(1:32,index2slag)) #table of all S_t1,... combinations
ref.s #structure of this reference table is crucial in step5!

slag2index<-function(slag){
  # convert a logical vector representing a into its integer form
  exponents<-rev(2^seq.int(0,length(slag)-1))
  #print(exponents)
  return(sum(exponents*slag))
}

slag2index(c(0,1,1)) #should return 3

lim.pi<-function(p,q){
  # Return the limiting value of the p,q combination chain S(t).
  return((1-q)/(1-p+1-q))
}

s.trans.prob<-function(snew,sold,p,q){
  # return the log-probability of a transition in S(t).
  # snew vector of S(t) values.
  # sold vector of corresponding S(t-1) values.
  log.probs<-rep(NA,length(snew))
  log.probs[snew==T & sold==T]<-log(p)
  log.probs[snew==F & sold==T]<-log(1-p)
  log.probs[snew==F & sold==F]<-log(q)
  log.probs[snew==T & sold==F]<-log(1-q)
  return(log.probs)
}

step0<-function(p,q){
  # return the starting probabilities at time 0
  probs<-matrix(NA,nrow=16,ncol=4)
  pival<-lim.pi(p,q)
  probs[ref.s[1:16,5]==T,4]<- log(pival)
  probs[ref.s[1:16,5]==F,4]<- log(1-pival)
  for(i in c(3,2,1)){
    probs[,i]<-s.trans.prob(ref.s[1:16,(i+1)],ref.s[1:16,(i+2)],p,q)
  }
  return(rowSums(probs))
}

iprobs<-step0(p,q)
sum(exp(iprobs)) #sanity check, should sum to one.

step1<-function(s.oldprobs,p,q){
  # s.oldprobs is the initializing 16-vector of log-probs from previous iteration
  # returns a 32-vector of log-probabilities (joint with s(t))
  transprobs<-s.trans.prob(ref.s[,1],ref.s[,2],p,q)
  return(transprobs+rep(s.oldprobs,2)) #joint probabilities of S(t),S(t-1)...
}

step2<-function(t,s.jointprobs,y,theta){
  # returns joint conditional density of y(t) and s.jointprobs (32-vector)
  # theta=(alpha1,alpha0,p,q,sigma,phi1,phi2,phi3,phi4)
  alpha1<-theta[1]
  alpha0<-theta[2]
  p<-theta[3]
  q<-theta[4]
  sigma<-theta[5]
  phi<-theta[6:9]
  deviations<-y[t]-alpha1*ref.s[,1]-alpha0 #32-vector
  means<-rep.int(0,32) #32-vector
  for(j in 1:4){
    if(t-j < 1) break #avoid exceeding the index boundaries
    means<-means+phi[j]*(y[t-j]-alpha1*ref.s[,(j+1)]-alpha0)
  }
  conditionals.y<-dnorm(deviations,mean=means,sd=sigma,log=TRUE)
  return(s.jointprobs+conditionals.y)
}

step3<-function(sy.jointprobs){
  # Given the 32-vector from step2, return the log total probability of y[t] (a scalar).
  sy.jointprobs<-exp(sy.jointprobs) #convert to probabilities to allow sum
  return(log(sum(sy.jointprobs)))
}

step4<-function(sy.jointprobs,y.totalprob){
  # first argument is from step2, second argument is from step3
  return(sy.jointprobs-y.totalprob)
}

step5<-function(s.fullprobs){
  #Takes the 32-vector output of step4, and consolidates it into a 16-vector by aggregating all the entries having the same 5th column indicator in the reference table of s(t). This is a way of marginalizing the contribution of the s[t-r] variable.
  # We have to be careful here about re-assigning to the correct (new) index when the frame shifts up by one. Fortunately the use of binary indices makes this possible (see the ref.s table).
  s.fullprobs<-exp(s.fullprobs)
  consol.probs<-rep(NA,16)
  for(i in 1:16){
    consol.probs[i]<-sum(s.fullprobs[(2*i-1):(2*i)])
  }
  return(log(consol.probs)) #can be reused in step1 for next t.
}

log.likelihood<-function(theta,y){
  #theta=(alpha1,alpha0,p,q,sigma,phi1,phi2,phi3,phi4).
  N<-length(y)
  p<-theta[3]
  q<-theta[4]
  s.oldprobs<-step0(p,q)
  logl<-0
  for(t in 1:N){
    s.jointprobs<-step1(s.oldprobs,p,q) 
    sy.jointprobs<-step2(t,s.jointprobs,y,theta)
    y.totalprob<-step3(sy.jointprobs)
    logl<-logl+y.totalprob
    s.fullprobs<-step4(sy.jointprobs,y.totalprob)
    s.oldprobs<-step5(s.fullprobs)
  }
  return(logl)
}