PYMC3 'FreeRV' object does not support item assignment erro

Question

I want to specify an array of random variables where random variable at index i has a sigma that is a function of random variable i-1, for i=2, 3, 4, ... N. I cannot see how to do this in pymc3.Here is my shot at it:

import pymc3 as pm

with pm.Model() as model:
    mu=pm.Normal('mu',-10,5)
    phi=pm.Beta('phi',1,1)
    tausq=pm.InverseGamma('tausq', 0,10000)
    tau=pm.Deterministic('tau', pm.math.sqrt(tausq))
    h0=pm.Normal('h0', 0,10000)
    h=pm.Normal('h',mu+phi*(h0-mu),tau, shape=len(returns))
    for i in range(2,len(returns)):
        h[i]=pm.Normal('h'+str(i),mu+phi*(h[i-1]-mu), tau)
    for i in range(1,len(returns)):
        r[i]=pm.Normal('r'+str(i),0, pm.math.exp(h[i]/2), observed=returns[i])
    trace_mdl=pm.sample(draws=2000,tune=2000,cores=10,init='adapt_diag')

This fails with error:

Traceback (most recent call last):

  File "<ipython-input-7-011e4fcf40bd>", line 11, in <module>
    h[i]=pm.Normal('h'+str(i),mu+phi*(h[i-1]-mu), tau)

TypeError: 'FreeRV' object does not support item assignment

I kinda understand that Pymc3 has the model specification component that I'm using here and perhaps I should be using the FreeRV class but that's as far as I can get. FWIW, this is how the original paper says its done in rstan:

parameters {
  real mu;
  real<lower=0,upper=1> phi;
  real<lower=0.01> tausq;
  real h0;
  vector[T] h;
}
transformed parameters {
  real<lower=0> tau;
  tau <- sqrt(tausq);
}

model {
  mu ~ normal(-10,5);
  phi ~ beta(1, 1);
  tausq ~ inv_gamma(2.5, 0.025);
  h0 ~ normal(0, 10000);
  h[1] ~ normal(mu + phi*(h0-mu), tau);
  for (t in 2:T)
    h[t] ~ normal(mu + phi*(h[t-1]-mu), tau);
  for (t in 1:T)
    y[t] ~ normal(0,exp(h[t]/2));
}