Derivation of the distribution of $\hat{\phi}=[\hat{\phi}_1, \cdots, \hat{\phi}_p]$ in AR(p) models Background
Consider the following AR($p$) model:
$$
\dot{X_t} = \phi_1 \dot X_{t-1} + \phi_2 \dot X_{t-2} + \cdots + \phi_p \dot X_{t-p} + \epsilon
$$
where $\dot{X} := X - \mu = X - \mathbb{E}(X)$, and $\epsilon \overset{iid}{\sim} N(0, \sigma_\epsilon^2) $.

Problem

Derive that the asymptotic distribution of $\hat{\phi}:=[\hat{\phi}_1, \cdots, \hat{\phi}_p]$ follows,
$$
\sqrt{n}(\hat{\phi} - \phi) \sim N(0, \sigma_\epsilon^2 \Gamma^{-1})
$$

where $\Gamma := \begin{bmatrix}
\gamma_0 & \gamma_1 & \gamma_2 & \cdots & \gamma_{k-1} \\
\gamma_1 & \gamma_0 & \gamma_1 & \cdots & \gamma_{k-2} \\
\gamma_2 & \gamma_1 & \gamma_0 & \cdots & \gamma_{k-3} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
\gamma_{k-1} & \gamma_{k-2} & \gamma_{k-3} & \cdots & \gamma_0\\
\end{bmatrix}$
with $\gamma_k := \text{Cov}(X_t, X_{t-k})$ : the autocovariance of lag $k$.

Try
To estimate $\phi_1, \cdots, \phi_p$, we usually consider the Yule-Walker equation,
$$
\mathbf{\rho} := \begin{bmatrix} \rho_1 \\ \rho_2 \\ \vdots \\ \rho_p \end{bmatrix}  = \begin{bmatrix}
1 & \rho_1 & \rho_2 & \cdots & \rho_{k-1} \\
\rho_1 & 1 & \rho_1 & \cdots & \rho_{k-2} \\
\rho_2 & \rho_1 & 1 & \cdots & \rho_{k-3} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
\rho_{k-1} & \rho_{k-2} & \rho_{k-3} & \cdots & 1 \\
\end{bmatrix} \begin{bmatrix} \phi_1 \\ \phi_2 \\ \vdots \\ \phi_p \end{bmatrix} = \mathbf{P}\mathbf{\phi}
$$
where $\phi_i$ : the coefficient of the model, and $\rho_k$ is the autocorrelation function(ACF) of lag $k$.

I think we can use the fact that we can estimate $\mathbf{\phi}$
$$
\hat{\mathbf{\phi}} = \hat{\mathbf{P}}^{-1}\hat{\mathbf{\rho}}
$$
where $\hat{\mathbf{\rho}}=[\hat{\rho_1}, \cdots, \hat{\rho_k}]^T$, which is the estimated version of $\rho$.

We can estimate $\rho_k$ which is a function of lag $k$ given data $\mathbf{x}$, as follows.
$$
\hat{\rho}_k := \frac{ \sum_{t=1}^{n-k} (x_t - \bar{x})(x_{t+k} - \bar{x}) }{\sum_{t=1}^n (x_t - \bar{x})^2}
$$
where $\bar{x} := \left(\sum_{t=1}^n x_t \right)/n$.
But I cannot proceed to show the result I want,
$$
\sqrt{n}(\hat{\phi} - \phi) \sim N(0, \sigma_\epsilon^2 \Gamma^{-1})
$$
 A: It doesn't have to be thru Yule-Walker. You can use the OLS approach, namely,
$$X_t = \pmb X_{t-1}^T \pmb \phi + n_t , \qquad t = p+1,p+2,\ldots$$
where $\pmb X_{t-1} = [X_{t-1} \ldots X_{t-p} ]$and $\pmb \phi = [\phi_1 \ldots \phi_p] $.
Stacking all such equations, we get
$$\pmb y = \pmb X \pmb \phi + \pmb n$$ 
where $\pmb y = [X_{p+1} \ldots X_N]^T$ and $$\pmb X = \begin{bmatrix} 
\pmb X_{p}^T\\
\pmb X_{p+1}^T\\
\vdots \\
\pmb X_{n-1}^T
\end{bmatrix} $$
Using OLS, we have 
$$\hat{\pmb{\phi}} = (\pmb X^T \pmb X)^{-1} \pmb X^T \pmb y = (\pmb X^T \pmb X)^{-1} \pmb X^T (\pmb X \pmb \phi + \pmb n) = \pmb \phi + (\pmb X^T \pmb X)^{-1} \pmb X^T \pmb n $$
Notice that from the above equation
$$\mathsf{E}( \hat{\pmb{\phi}} ) = \pmb \phi +(\pmb X^T \pmb X)^{-1} \underbrace{\mathsf{E} \pmb X^T \pmb n}_0 \pmb \phi$$.
and
$$\mathsf{var}(\hat{\pmb{\phi}}) = \mathsf{E}( (\hat{\pmb{\phi}} - \pmb\phi)(\hat{\pmb{\phi}} - \pmb\phi)^T ) =(\pmb X^T \pmb X)^{-1} \pmb X^T \mathsf{E} (\pmb n \pmb n^T) \pmb X (\pmb X^T \pmb X)^{-1} $$
assuming white noise with same variance across time $\sigma_\epsilon^2$, we get
$$\mathsf{var}(\hat{\pmb{\phi}}) = \mathsf{E}( (\hat{\pmb{\phi}} - \pmb\phi)(\hat{\pmb{\phi}} - \pmb\phi)^T ) =(\pmb X^T \pmb X)^{-1} \pmb X^T \sigma_\epsilon^2 \pmb I \pmb X (\pmb X^T \pmb X)^{-1} = \sigma_\epsilon^2 (\pmb X^T \pmb X)^{-1}(\pmb X^T \pmb X)(\pmb X^T \pmb X)^{-1}=\sigma_\epsilon^2 (\pmb X^T \pmb X)^{-1}$$
Now assuming your AR(p) is strictly stationary and ergodic and that $\mathsf{E}(X_t^4)$ exists, then the central limit theorem applies
$$\sqrt{n}( \hat{\pmb{\phi}} - \pmb{\phi}) \longrightarrow N(0,\mathsf{var}(\hat{\pmb{\phi}}))$$ 
Note that $(\pmb X^T \pmb X)^{-1} = \Gamma^{-1}$ 
