How to estimate CAViaR (Engle and Manganelli 2004) using non linear quantile regression?

I am trying to replicate results from Engle and Manganelli (2004). The following is one of their specifications, $q_t(\theta)=\gamma_0+\gamma_1q_{t-1}(\theta)+\alpha|r_{t-1}|$, $q$ is the quantile of return distribution, $r$ is the return.

I do not know how to use quantile regression to estimate this process, since we have quantiles on both sides of equation. Any suggesions?

You have to minimize RQ = - ( $$I$$($$r$$ < - $$q$$ ) - $$\theta$$ ) * ( $$r$$ + $$q$$ ). Where $$r$$,$$q$$ are the vectors of returns and quantiles found (for a given $$\gamma$$0, $$\gamma$$1, $$\alpha$$ ), * is the matrix multiplication, and $$I$$( ) is the indicator function. Note that because we have an autoregressive function, you must have a first $$q$$ (VaR) which in the paper is the choosen $$\theta$$ quantile of the first 300 observations. The minimization problem is solved by finding which coefficents ($$\gamma$$0,$$\gamma$$1 and $$\alpha$$ ) is the ones that minimizes RQ. In the paper this is done by generating 10000 random coefficients combinations, then using the top 10 $$a$$ results in a optimizer ( he uses 'fminsearch' and 'fminunc' functions from MATLAB, in R you can use optim() ) and reiterating 5 times to refine the results.
$$a$$ - for the simetric absolute CAViaR