Simulate a model using the specified decision rule.

    Parameters
    ----------

    model: NumericModel

    dr: decision rule

    s0: ndarray
        initial state where all simulations start

    driving_process: ndarray
        realization of exogenous driving process (drawn randomly if None)

    N: int
        number of simulations
    T: int
        horizon for the simulations
    seed: int
        used to initialize the random number generator. Use it to replicate
        exact same results among simulations
    discard: boolean (False)
        if True, then all simulations containing at least one non finite value
        are discarded

    Returns
    -------
    xarray.DataArray:
        returns a ``T x N x n_v`` array where ``n_v``
           is the number of variables.

def simulate(model, dr, N=1, T=40, s0=None, i0=None, m0=None,
        driving_process=None, seed=42, stochastic=True):
    '''
    Simulate a model using the specified decision rule.

    Parameters
    ----------

    model: NumericModel

    dr: decision rule

    s0: ndarray
        initial state where all simulations start

    driving_process: ndarray
        realization of exogenous driving process (drawn randomly if None)

    N: int
        number of simulations
    T: int
        horizon for the simulations
    seed: int
        used to initialize the random number generator. Use it to replicate
        exact same results among simulations
    discard: boolean (False)
        if True, then all simulations containing at least one non finite value
        are discarded

    Returns
    -------
    xarray.DataArray:
        returns a ``T x N x n_v`` array where ``n_v``
           is the number of variables.
    '''

    if isinstance(dr, AlgoResult):
        dr = dr.dr

    calib = model.calibration
    parms = numpy.array(calib['parameters'])

    if s0 is None:
        s0 = calib['states']

    n_x = len(model.symbols["controls"])
    n_s = len(model.symbols["states"])

    s_simul = numpy.zeros((T, N, n_s))
    x_simul = numpy.zeros((T, N, n_x))

    s_simul[0, :, :] = s0[None, :]

    # are we simulating a markov chain or a continuous process ?
    if driving_process is not None:
        m_simul = driving_process
        sim_type = 'continuous'
        m0 = m_simul[0,:,:]
        x_simul[0,:,:] = dr.eval_ms(m0, s0)
    else:
        if isinstance(dr.exo_grid, UnstructuredGrid):
            if i0 is None:
                i0 = 0
            dp = model.exogenous.discretize() ## TODO (nothing guarantee the processes match)
            m_simul = dp.simulate(N, T, i0=i0, stochastic=stochastic)
            i_simul = find_index(m_simul, dp.values)
            sim_type = 'discrete'
            m0 = dp.node(i0)
            x0 = dr.eval_is(i0, s0)

        else:
            process = model.exogenous ## TODO (nothing guarantee the processes match)
            m_simul = process.simulate(N, T, m0=m0, stochastic=stochastic)
            sim_type = 'continuous'
            if m0 is None:
                m0 = model.calibration["exogenous"]
            x0 = dr.eval_ms(m0, s0)
            x_simul[0, :, :] = x0[None, :]

    fun = model.functions
    f = model.functions['arbitrage']
    g = model.functions['transition']

    numpy.random.seed(seed)
    from dolo.misc.dprint import dprint
    mp = m0
    for i in range(T):
        m = m_simul[i,:,:]
        s = s_simul[i,:,:]
        if sim_type=='discrete':
            i_m = i_simul[i,:]
            xx = [dr.eval_is(i_m[ii], s[ii,:]) for ii in range(s.shape[0])]
            x = np.row_stack(xx)
        else:
            x = dr.eval_ms(m, s)

        x_simul[i,:,:] = x
        ss = g(mp, s, x, m, parms)
        if i < T-1:
            s_simul[i + 1, :, :] = ss
        mp = m

    if 'auxiliary' not in fun:  # TODO: find a better test than this
        l = [s_simul, x_simul]
        varnames = model.symbols['states'] + model.symbols['controls']
    else:
        aux = fun['auxiliary']
        a_simul = aux(
            m_simul.reshape((N * T, -1)),
            s_simul.reshape((N * T, -1)),
            x_simul.reshape((N * T, -1)), parms)
        a_simul = a_simul.reshape(T, N, -1)
        l = [m_simul, s_simul, x_simul, a_simul]
        varnames = model.symbols['exogenous'] + model.symbols['states'] + model.symbols[
            'controls'] + model.symbols['auxiliaries']

    simul = numpy.concatenate(l, axis=2)

    import xarray as xr
    data = xr.DataArray(
            simul,
            dims=['T','N','V'],
            coords={'T': range(T), 'N': range(N), 'V': varnames}
        )

    return data