Proof expression for the autocovariance function of AR(1) The representation for the model AR(1) is the following:
$Y_t=c+ϕY_{t-1}+ε_t$
where $c=(1-ϕ)μ$ ($c$ is a constant).

I want to understand the calculations that there are behind the general formula of the autocovariance of AR(1), which is $γ(h)=\operatorname{Var}(Y_t )⋅ϕ^{|h|} $

So far, I did the following steps - I started with $γ(1)$:
$\operatorname{Cov}(Y_t,Y_{t-1})=γ(1)=$
$=\operatorname{Cov}(ϕY_{t-1}+ε_t, ϕY_{t-2}+ε_t)=$
$=\operatorname{Cov}(ϕY_{t-1}, ϕY_{t-2}) + \operatorname{Cov}(ϕY_{t-1}, ε_t) + \operatorname{Cov}(ε_t, ϕY_{t-2}) +  \operatorname{Cov}(ε_t, ε_t)
$

$γ(1)=ϕ^2γ(1) + ???+???+σ^2$

As you can see, from this point I can't continue because I don't know which are the values of $\operatorname{Cov}(ϕY_{t-1}, ε_t)$ and $\operatorname{Cov}(ε_t, ϕY_{t-2})$

Any assistance will be much appreciated. Thank you in advance.
 A: Starting from what you have provided:
$y_{t} = c + \phi y_{t-1} + \epsilon_{t} \tag{1}$ 
Where $c = (1 - \phi) \mu$ 

We can rewrite $(1)$ as:   
\begin{array}
\ y_{t} & = & c + \phi y_{t-1} + \epsilon_{t} \\
& = & (1 - \phi) \mu + \phi y_{t-1} + \epsilon_{t} \\
& = & \mu - \phi \mu + \phi y_{t-1} + \epsilon_{t} \\
\end{array}
Then,   
$y_{t} - \mu = \phi (y_{t-1} - \mu) + \epsilon_{t} \tag{2}$ 
If we let $\tilde{y_{t}} = y_{t} - \mu$, then equation $(2)$ can be writen as:   
$\tilde{y}_{t} = \phi \tilde{y}_{t-1} + \epsilon_{t} \tag{3}$ 

Variance 
The variance of $(3)$ is obtained by squaring the expression and taking expectations, which ends in:    
\begin{array}
\ \tilde{y}_{t}^2 & = & (\phi \tilde{y}_{t-1} + \epsilon_{t})^2 \\
& = & (\phi \tilde{y}_{t-1})^2 + 2 \phi \tilde{y}_{t-1} \epsilon_{t} + (\epsilon_{t})^2 \\
& = & \phi^{2} \tilde{y}_{t-1}^{2} + 2 \phi \tilde{y}_{t-1} \epsilon_{t} + \epsilon_{t}^2 
\end{array} 
Now take the expectation:   
$E(\tilde{y}_{t}^2) = \phi^{2} E(\tilde{y}_{t-1}^{2}) + 2 \phi E(\tilde{y}_{t-1} \epsilon_{t}) + E(\epsilon_{t}^2)$ 
Here we will call:   


*

*$\sigma_{y}^{2}$ is the variance of the stationary process.   

*The second term in the right-hand side of the equation is zero because $\tilde{y}_{t-1}$ and $\epsilon_{t}$ are independent and both have null expectation.    

*The last term in the right is the variance of the innovation, denoted as $\sigma^{2}$ (note that there is no subscript for this).    


Finally, 
$\sigma_{y}^{2} = \phi^{2} \sigma_{y}^{2} + \sigma^{2}$ 
If we solve for the variance of the process, namely $\sigma_{y}^{2}$, we have:   
$\sigma_{y}^{2} = \frac{\sigma^2}{1 - \phi^2} \tag{4}$ 

Autocovariance 
We are going to use the same trick we use for formula $(3)$. The autocovariance between observations separated by $h$ periods is then:  
\begin{array}
\ \gamma_{h} & = & E [(y_{t - h} - \mu) (y_{t} - \mu)] \\
& = & E[(\tilde{y}_{t - h}) (\tilde{y}_{t})] \\
& = & E[\tilde{y}_{t - h} (\phi \tilde{y}_{t - 1} + \epsilon_{t}) \\
\end{array} 
The innovations are uncorrelated with the past values of the series, then $E[\tilde{y}_{t-h} \epsilon_{t}] = 0$ and we are left with:  
$\gamma_{h} = \phi \gamma_{h-1} \tag{5}$ 
For $h = 1, 2, \ldots$  and with $\gamma_{0} = \sigma_{y}^2$ 

For the particular case of the $AR(1)$, equation $(5)$ becomes:  
$\gamma_{1} = \phi \gamma_{0}$ 
And using the result from equation $(4)$: $\gamma_{0} = \sigma_{y}^{2} = \frac{\sigma^2}{1 - \phi^2}$ we end up with    
$\gamma_{1} = \frac{\sigma^2}{1 - \phi^2} \phi$

Original source: Andrés M. Alonso & Carolina García-Martos slides. Available here: http://www.est.uc3m.es/amalonso/esp/TSAtema4petten.pdf 
A: Let's write $\gamma(1)$:
$$\begin{align}\gamma(1) &= cov(Y_t,Y_{t-1})=cov(c+\phi Y_{t-1}+\epsilon_t,Y_{t-1})\\&=cov(c,Y_{t-1})+\phi cov(Y_{t-1},Y_{t-1})+cov(\epsilon_t,Y_{t-1}) \\&=\phi \gamma(0)\end{align}$$ 
since $cov(c,Y_{t-1})=cov(\epsilon_t,Y_{t-1})=0$, (i.e. past output is independent from future input). 
Similarly, $\gamma(2)=cov(Y_t,Y_{t-2})=cov(\phi Y_{t-1}+c+\epsilon_t,Y_{t-2})=\phi \gamma(1)=\phi^2\gamma(0)$. 
If we continue this way, we get $cov(Y_t,Y_{t-h})=\phi \gamma(h-1)=...\phi^h\gamma(0)$, where $h\geq0$. Generalizing for negative $h$ yields $\gamma(h)=\phi^{|h|}\gamma(0)$, where $\gamma(0)=var(Y_t)$.
P.S. all of this analysis assumes $\epsilon_t$ is WSS, therefore $y_t$ from LTI filtering property.
