Consider the $m$-dimensional VAR process
$${\bf x}_t = \sum_{l=1}^{P} A_l{\bf x}_{t-l} + {\bf e}_t$$
where the componenets of ${\bf e}_t$ are spatially and temporally independent and follow a generalised Gaussian distribution i.e. $e_{i,t} \sim GGD(\alpha_i, \beta_i)$.
Given my data $X$, I am trying to estimate the VAR coefficients $A_l$ and GGD parameters $\alpha_i, \beta_i$. I have tried to use Newton-Raphson optimisation with the vector
$$[\alpha_1, ..., \alpha_M, \beta_1, ..., \beta_M, A_1(1,1), ..., A_1(1,M), A_1(2,1), ..., A_P(M,M)]^T$$
but it is painfully slow and not very robust.
What would be the best way to tackle this model fitting problem?
