Polinomio característico de un VAR

[Víctor de Buen Remiro]

Se solicita incoporar a VarModel el cálculo del polinomio característico de un modelo VAR y el cálculo de sus raíces. Esto es necesario para diagnosticar la estabilidad del modelo VAR. Se aporta el código que lo implementa:

Matrix param = varEst::_.B;
Real n = Rows(param);
Set PHI = 
{
  Set aux = For(1,pmax,Set(Real k)
  {
    Matrix phi = Sub(param,1,(k-1)*n+1,n,n);
    For(1,n,Set(Real i)
    {
      For(1,n,Polyn(Real j)
      {
        MatDat(phi,i,j)*B^k
      })
    })
  });
  For(1,n,Set(Real i)
  {
    For(1,n,Polyn(Real j)
    {
      If(i==j,1,0)-
      SetSum(For(1,pmax,Polyn(Real k)
      {
        aux[k][i][j]
      })) 
    })
  })
};

Polyn PolynMatrixDeterminant(Set polMat)
{
  Real m = Card(polMat);
  WriteLn("TRACE [PolynMatrixDeterminant] polMat=\n"<<polMat);
  WriteLn("TRACE [PolynMatrixDeterminant] m="<<m);
  Case(
    m==1, polMat[1][1],
    m==2, polMat[1][1]*polMat[2][2]-polMat[1][2]*polMat[2][1], /*
    m==-3, 
    {
      +polMat[1][1]*polMat[2][2]*polMat[3][3]
      +polMat[1][2]*polMat[2][3]*polMat[3][1]
      +polMat[1][3]*polMat[2][1]*polMat[3][2]
      -polMat[1][1]*polMat[2][3]*polMat[3][2]
      -polMat[1][2]*polMat[2][1]*polMat[3][3]
      -polMat[1][3]*polMat[2][2]*polMat[3][1]
    }, */
    1==1,
    {
      SetSum(For(1,m,Polyn(Real k)
      {
        Polyn c = polMat[k][1]*(-1)^(k+1);
        Set rows = ExtractByIndex(polMat,Range(1,m,1)-[[k]]);
        Set adjunt = EvalSet(rows,Set(Set row)
        {
          For(2,m,Polyn(Real j) { row[j] })
        });
        c*PolynMatrixDeterminant(adjunt)
      }))
    }
  )
};
  

Polyn PHI.det = PolynMatrixDeterminant(PHI);

Set PHI.det.roots = 
{
  Matrix ar.roots = gsl_poly_complex_solve(PHI.det);
  For(1,Rows(ar.roots),Set(Real k)
  {
    Set aux = [[
      Real real = MatDat(ar.roots,k,1);
      Real imaginary = MatDat(ar.roots,k,2);
      Real module = Sqrt(real^2+imaginary^2)
    ]];
    Eval("Set PHI.det.Root_"<<k+"=aux")
  })
};