/* Figure 13.4 in POE */ 

clear
set obs 30                            
gen  time = _n
tsset time

mata
delta_10 = 0
delta_11 = 0.7
delta_12 = 0.2
delta_20 = 0
delta_21 = 0.3
delta_22 = 0.6

result = J(30,4,0)
result[1,1]=1
result[1,2]=0
result[1,3]=0 
result[1,4]=2


i = 2
while (i<=30) {
   result[i,1]= delta_10 + delta_11*result[i-1,1] + delta_12*result[i-1,2]
   result[i,2]= delta_20 + delta_21*result[i-1,1] + delta_22*result[i-1,2]
   result[i,3]= delta_10 + delta_11*result[i-1,3] + delta_12*result[i-1,4] 
   result[i,4]= delta_20 + delta_21*result[i-1,3] + delta_22*result[i-1,4]
   i++
}

result
  (void) st_addvar("double","yresptoy")
  (void) st_addvar("double","xresptoy")
  (void) st_addvar("double","yresptox")
  (void) st_addvar("double","xresptox")

  st_store(.,"yresptoy",result[.,1])
  st_store(.,"xresptoy",result[.,2])
  st_store(.,"yresptox",result[.,3])
  st_store(.,"xresptox",result[.,4])

  end

list

tsline yresptoy, title(response of y to y) saving(figure13pt4a, replace)
tsline xresptoy, title(response of x to y) saving(figure13pt4c, replace)
tsline yresptox, title(response of y to x) saving(figure13pt4b, replace)
tsline xresptox, title(response of x to x) saving(figure13pt4d, replace)
graph combine figure13pt4a.gph figure13pt4b.gph figure13pt4c.gph figure13pt4d.gph, saving(figure13.4, replace)
save figure13.4.dta, replace