Let’s take a simple example to answer a broad but interesting question:
Imagine that we have a daily return serie denoted $r_{t}$ ( which is assumed to be stationary) and let's take a little time to define main concepts :
Mean Process (First moment process)
The unconditional mean of $r_{t}$ denoted $u$ is just its expectation $E(r_{t})$. It is not time varying. You can compute it directly using the expectation formula.
The conditional mean process refers to the expectation of the serie at time $t$ given previous information: $E_{t}(y_{t}|\Omega_{t-1})$. It is time varying and that is the reason way we write it using a time subscript: $u_{t}$. This process is usually estimated using autoregressive–moving-average (ARMA) models. The intuition is that we can detect some autocorrelations in the returns series (ex: if day one is up, the day after has more probability to be down...it is an example)
So far so good, we may assume that we can compute a simple "average" return (unconditional mean) or a time-varying (i.e conditional) "average" return.
However usually people are also concerned with risk. If you know that your returns in average follow a process, you are likely also interested in uncertainty/risk. In finance, risk is usually approximated using the second moment (ie the variance).
Now let's jump to the variance part:
Conditional variance Process (Second moment)
Similarly that for the mean process, we are able to estimate the unconditional variance of our return serie using a simple variance formula $\sigma^{2} = Var(r_{t}) $.
Now imagine our return series exhibits "large changes followed by large changes..." during few days and goes back to its original unconditional variance level. We may realize that the variance is in fact time varying : we observe some "volatility clustering". In a same way that for the conditional mean process we can build a conditional variance process. To this end we use different tools : the Garch family models which allows us to model a time-varying variance : $\sigma_{t}^{2} = Var_{t}(r_{t}|\Omega_{t-1}) $. (Others models exist such as Stochastic volatility models).
Now we have defined the main concepts we can jump to your question :
How to calculate the conditional variance of a time series?
Intuition
Firstly we model the conditional mean process (using a ARMA,ARFIMA...) and subtract it from the original returns series to obtain the "return residuals" : $r_{t}-\mu_{t} =\epsilon_{t} = \sigma_{t} z_{t}$ where $z_{t}$ is an i.i.d process with $E_{t}(z_{t}) = 0$ and $Var_{t}(z_{t}) =1$. Note that the conditional variance of $\epsilon_{t}$ is equal to $\sigma_{t}^{2}$.
However since we know that the variance is time varying we also know that $\sigma_{t}^{2}$ has a time dependent structure and exhibits autocorrelations (so do the squares returns residuals). We can model it using GARCH class of models which can (very roughly) be seen as ARMA models for the conditional variance process.
Example of a Garch(1,1) : $\sigma_{t}^{2} = a + \alpha \epsilon_{t-1}^{2} + \beta \sigma_{t-1}^{2} $
Once we fit our conditional variance models we will be left with the conditional variance process $\sigma_{t}^{2} $.At this point we know the conditional variance process $\sigma_{t}^{2} $ and $\epsilon_{t}^{2}$. This allow us to obtain the final standardized residuals series $z_{t}$ which is i.i.d and equal to $\epsilon_{t}/\sigma_{t} = z_{t}$.
Estimation
How do we estimate it ?
The simplest way is to rely on the Maximum Likelihood Estimation (MLE) method. We need to assume a distribution for the $z_{t}$ (the final residuals). Since we know that these residuals are i.i.d it is easy to compute the log-likelihood for a given $z_{t}$ serie. (to be more precise the typical arguments for the likelihood function are $\epsilon_{t}^{2}$ and $\sigma_{t}^{2}$)
Example : If we assume a normal distribution for $z_{t}$ the log likelihood (assuming no constant) is given by : $LogLik = -\frac{1}{2} \sum_{1}^{T} \left[ \log(2 \pi) + \log(\sigma_{t}^{2}) + z_{t}^{2} \right] \qquad (= -\frac{1}{2} \sum_{1}^{T} \left[ \log(2 \pi) + \log(\sigma_{t}^{2}) + \frac{\epsilon_{t}^{2}}{\sigma_{t}^{2}} \right])$
But how can we practically obtain $z_{t}$ ? A solution is to use what we called "filters" takings as input the returns series and, based on a particular specifications (ex: arma(1,1)-garch(1,1)), returning $\sigma_{t}^{2}$. By "filtering" we mean that we applied the autoregressive framework (recursive algorithm on both the mean and variance) on the input serie (the return serie) to obtain as output the $ z_{t}$ .
An example ? : see this very nice post (see Solution added by the author) : Algorithm to fit AR(1)/GARCH(1,1) model of log-returns
Next we can use some maximization algorithms to find the parameters producing a $z_{t}$ serie which maximize the likelihood.
Ex: we run the filter with AR parameter = 0.1 ; next we try another value, and so on with all parameters to obtain the final parameters maximizing the likelihood.
Finally to obtain the standards errors of the estimated parameters we may use the Hessian.
Ok, but in practice ?
You can use "click to run" softwares to estimate (and much more) the conditional variance/mean processes. Matlab, R , and Ox (among others) have packages devoted to this estimation.
Example of packages :
Remarks
-This is a simplification example: models currently used in the literature are much more advanced, for instance the Arch-in-mean class of models add the conditional variance as an explanatory variable in the conditional mean process.
-You are not forced to use filters if you can compute directly the likelihood based on the parameters.
-In reality, the estimation part is far more difficult to do, as illustration the choice of starting values is a tricky part.
-If you want to recommend another software/package just add it in the comments.