#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

at = c(-0.3,0.6,0.9,0.2,0.1,-0.6,1.7,-0.9,-1.3,-0.6,-0.4,0.9,0.0,-1.4,-0.6) # a_0, a_1, to a_{14} 
Na = length(at)
Nam1 = Na-1

zt0 = c(19,20) # z_{-1}, z_0

# First compute the time series values using the difference equation form:
# 

# Part (a):
#
zt_a = zt0
for( ii in 1:Nam1 ){ # the index of z_t 
  newZ = zt_a[length(zt_a)] + at[ii+1] - 0.5 * at[ii] 
  zt_a = c( zt_a, newZ )
}
   
# Part (b):
#
zt_b = zt0
for( ii in 1:Nam1 ){ # the index of z_t 
  newZ = zt_b[length(zt_b)] + at[ii+1] - 0.2 * at[ii] 
  zt_b = c( zt_b, newZ )
}

# Part (c):
#
zt_c = zt0
for( ii in 1:Nam1 ){ # the index of z_t 
  newZ = 1.5*zt_c[length(zt_c)] -0.5*zt_c[length(zt_c)-1] + at[ii+1] 
  zt_c = c( zt_c, newZ )
}

# Part (d):
#
zt_d = zt0
for( ii in 1:Nam1 ){ # the index of z_t 
  newZ = 1.2*zt_d[length(zt_d)] -0.2*zt_d[length(zt_d)-1] + at[ii+1] 
  zt_d = c( zt_d, newZ )
}

# Part (e):
#
zt_e = zt0
for( ii in 1:Nam1 ){ # the index of z_t 
  newZ = 1.2*zt_e[length(zt_e)] -0.2*zt_e[length(zt_e)-1] + at[ii+1] - 0.5 * at[ii] 
  zt_e = c( zt_e, newZ )
}

# The indexes of t in the representations of a_t and z_t above:
a_t_index =  0 : 14 
z_t_index = -1 : 14 

# Now compute the time series values using the recursive forms:
#

# Part (a):
#
for( ii in 12:14 ){ # we want to compute z_{ii} 
  Nz = ii + 2 - 1 # add 2 to use z_{-1} and z_0 index (subtract one to look at all z_{t} with t < ii)
  pij = rev( 0.5**( 1:Nz ) ) 
  zNew = at[ii+1] + sum( zt_a[1:Nz] * pij ) # add 1 from at to remove 0 index
  print(sprintf("(a): index t= %4d: difference equation= %6.4f recursive= %6.4f",ii,zt_a[ii+2], zNew))
  #zt_a[ii+2] = zNew # whether we want to use our most recent estimate in the next iteration
}

# Part (b):
#
for( ii in 12:14 ){ # we want to compute z_{ii} 
  Nz = ii + 2 - 1 # add 2 to use z_{-1} and z_0 index (subtract one to look at all z_{t} with t < ii)
  pij = 4*rev( 0.2**( 1:Nz ) ) 
  zNew = at[ii+1] + sum( zt_b[1:Nz] * pij ) # add 1 from at to remove 0 index
  print(sprintf("(b): index t= %4d: difference equation= %6.4f recursive= %6.4f",ii,zt_b[ii+2], zNew))
  #zt_b[ii+2] = zNew # whether we want to use our most recent estimate in the next iteration
}

# Part (c):
#
for( ii in 12:14 ){ # we want to compute z_{ii} 
  Nz = ii + 2 - 1 # add 2 to use z_{-1} and z_0 index (subtract one to look at all z_{t} with t < ii)
  zm1 = zt_c[Nz]
  zm2 = zt_c[Nz-1] 
  zNew = at[ii+1] + 1.5*zm1 - 0.5*zm2 # add 1 from at to remove 0 index
  print(sprintf("(c): index t= %4d: difference equation= %6.4f recursive= %6.4f",ii,zt_c[ii+2], zNew))
  #zt_c[ii+2] = zNew # whether we want to use our most recent estimate in the next iteration
}

# Part (d):
#
for( ii in 12:14 ){ # we want to compute z_{ii} 
  Nz = ii + 2 - 1 # add 2 to use z_{-1} and z_0 index (subtract one to look at all z_{t} with t < ii)
  zm1 = zt_d[Nz]
  zm2 = zt_d[Nz-1] 
  zNew = at[ii+1] + (6./5)*zm1 - (1./5)*zm2 # add 1 from at to remove 0 index
  print(sprintf("(d): index t= %4d: difference equation= %6.4f recursive= %6.4f",ii,zt_d[ii+2], zNew))
  #zt_d[ii+2] = zNew # whether we want to use our most recent estimate in the next iteration
}

# Part (e):
#
theta = 0.5
phi = 0.2 
for( ii in 12:14 ){ # we want to compute z_{ii} 
  Nz = ii + 2 - 1 # add 2 to use z_{-1} and z_0 index (subtract one to look at all z_{t} with t < ii)
  pi1 = phi + ( 1 - theta )
  pij = c( rev( (theta - phi)*(1-theta)*theta**( 0:(Nz-2) ) ), pi1 ) 
  zNew = at[ii+1] + sum( zt_e[1:Nz] * pij ) # add 1 from at to remove 0 index
  print(sprintf("(e): index t= %4d: difference equation= %6.4f recursive= %6.4f",ii,zt_e[ii+2], zNew))
  #zt_e[ii+2] = zNew # whether we want to use our most recent estimate in the next iteration
}