Newcomer question: How does the GARCH recursive formula actually work? So, I have a some experience with standard econometrics, and I also understand GARCH's basic concept, but I can't figure out what actually goes into its model.
So we have the standard GARCH(1,1) recursive equation:
$$\sigma_{t}^{2} = \alpha_0 + \alpha (\sigma_{t-1}^{2}) + \beta (\sigma_{t-1}^{2}\eta_{t-1}^2) $$
If I have returns of something, let's say the monthly S&P 500 for 10 years, what do I actually input into this equation? I have read here that your input needs to be a de-trended and de-seasonalized dataset; ok, so we apply some time-series regression to the returns first, which makes sense. What then? Do we apply the output of the time-series, or just the time-series' residuals to the equation, and where? Into the $\sigma$? But that isn't $\sigma$ at all, but a time-series.
I'm just very confused on how the function is supposed to work, guys, and would appreciate a step by step introduction from having one's data, to getting a GARCH model.
 A: In order to clarify some sources of potential confusion, I will use a slightly different notation, which will hopefully make it easier to follow through my answer. Also, I think @Johan's answer is very good, but I just wanted to add a few more pedagogical aspects, hence my attempt.
Consider the return series $r_t$ of your chosen asset. Typically, in an econometric framework, this is modeled as $r_t = \mu_t + \varepsilon_t$, where $\mu_t$ is the trend component, and $\varepsilon_t$ is the noise. A not so unreasonable assumption here (which I am doing only for simplification and does not change the result dramatically), as @Johan mentions, is to consider $\mu_t = \mu$, in which case it will very likely be statistically identical to $0$ (for reasons as to why this is assumed, check out the efficient markets hypothesis).
This assumption simplifies our model to $r_t = \varepsilon_t$, which we take one step further by expressing it as $\varepsilon_t = \sigma_t Z_t$, where $\sigma_t^2$ is the conditional variance of $r_t$, i.e. $\sigma_t^2 = \mathbb{V}[r_t\mid \mathcal{F}_{t-1}]$, and $Z_t$ are iid with $\mathbb{E}[Z_t] = 0$ and $\mathbb{V}[Z_t] = 1$.
Now, the topic of concern here is the way $\sigma_t^2$ evolves over time. In the GARCH(1,1) setting, this is described by
$$ \sigma_t^2 = \omega + \alpha \varepsilon_{t-1}^2 + \beta\sigma_{t-1}^2 = \omega + \alpha(\sigma_{t-1} Z_{t-1})^2 + \beta\sigma_{t-1}^2 $$
where $\omega>0$ the constant term, $\alpha$ and $\beta$ the multipliers of the ARCH & GARCH terms respectively (notice here that I have switched the order of these parameters in comparison to your original post).
To quote the "Quantitative Risk Management: Concepts, Techniques, and Tools" book by McNeil, Frey, and Embrechts (which I suggest you read in case you are not satisfied with the answers presented here), $\sigma_t$ is recursively defined in terms of $\sigma_{t-1}$, which one typically initializes by using the unconditional variance of the entire sample, i.e. $\sigma_0 = \mathbb{V}[r]$.
In the MATLAB code below, you can see how you can create the log-likeligood $f$ for a GARCH(1,1) model based on the data $r$, with $Z_t\stackrel{iid}{\sim} \mathcal{N}(0,1)$, and the small change in notation $h_t = \sigma_t^2$:
function f = garch(a, r)
T = length(r);
omega = a(1); % i.e. the constant term
a = a(2); % i.e. the multiplier of the ARCH term
b = a(3); % i.e. the multiplier of the GARCH term

h = [var(r); zeros(T-1,1)]; % vector of conditional variances
Z = [r(1)/sqrt(h(1)); zeros(T-1,1)]; % vector of innovations
v = [log(2*pi) + log(h(1)) + (Z(1))^2; zeros(T-1,1)]; % vector of individual log-likelihood contributions

for t = 2:T
    h(t) = omega + a*(r(t-1))^2 + b*(h(t-1))^2;
    Z(t) = r(t)/sqrt(h(t));
    v(t) = log(2*pi) + log(h(t)) + (Z(t))^2;
end
f = -0.5*sum(v);
end

You would then pass this function as the target function to be optimized w.r.t the vector a of parameters, and whose input would be the data r.
Please note I had to modify this function from a different GARCH model I had worked on in the past, and I don't have MATLAB currently installed so I can't test it, but hopefully it conveys the point I was trying to make.
EDIT: the assumption that $\mu_t = 0$ is not necessarily correct (it could even be flat out wrong, depending on the data), but was made here to simplify the operations. If the mean equation is not 0 (or if we don't want it to be), then the difference would simply be going from $Z_t = r_t/\sigma_t$ to $Z_t = (r_t - \mu_t)/\sigma_t$, and of course, any additional parameters in the expression for $\mu_t$ would need to be included accordingly in the log-likelihood function for the maximum likelihood estimation of the whole model's parameters.
A: Assume that we want to model the volatility of SP500 returns, $r_t$. 
A GARCH type model consists of both a mean and GARCH equation. The mean equation can be defined as
$$ r_t = \mu_t + \sigma_t \eta_t$$ 
where $\eta_t$ is iid with zero mean and unit variance and $\mu_t$ represents the mean dynamics. Often just assumed to be constant, but could be more advanced depending on the data. 
The idea is that the volatility of the error term in the mean equation is time-varying modelled with the GARCH equation. 
Lets assume that $\mu_t = \mu$, then we put demeand returns into the GARCH equation
$$\sigma_t^2 = \alpha_0 + \alpha_1\sigma_{t-1}^2+\beta (r_{t-1}-\mu)^2$$
We can think about these demeand returns as the signal about volatility. We do not know $\sigma_0^2$, why we need to initialize it (e.g. with $(r_0-\mu)^2$ or the unconditional variance).
