function [mcode,MaxAbsError,BMatrix,SMatrix,ModelEigenvalues] = AiMSolver(ModelParameters,NumEq,NumLag,NumLead,AIMTolerance)
% AiMSolver: Attempts to solve the model using AiM
%
% USAGE:
%
%       [mcode,MaxAbsError,BMatrix,SMatrix,ModelEigenvalues] = AiMSolver(ModelParameters,NumEq,NumLag,NumLead,AIMTolerance)
%
%
% REQUIRED INPUTS:  ModelParameters (structure) contains numeric values for the model parameters
%
%                   NumEq (integer), the number of equations of the model.
%
%                   NumLag (integer), the number of lags.
%
%                   NumLead (integer), the number of leads.
%
%                   AIMTolerance (positive integer) that determines the tolerance for the AiM solver.
%
% REQUIRED OUTPUTS: mcode (integer). Indicates whether the solution of the model is unique or not.
%                         mcode is 1 when unique, 0 if not computed, and various other values when
%                         there are problems with the solution. Specifically:
%
%                              mcode = 2:      Roots not correctly computed by real_schur
%                              mcode = 3:      Too many big roots
%                              mcode = 35:     Too many big roots, and q(:,right) is singular
%                              mcode = 4:      Too few big roots
%                              mcode = 45:     Too few big roots, and q(:,right) is singular
%                              mcode = 5:      q(:,right) is singular
%                              mcode = 61:     Too many exact shiftrights
%                              mcode = 62:     Too many numeric shiftrights
%                              mcode = 7:      The a matrix has infinite or NaN values. Eigenvalues
%                                              cannot be calculated.
%                              mcode = 8:      The function "compute_aim_matrices" returns complex
%                                              numbers.
%                              else            Return code not properly specified
%
%                   MaxAbsError (non-negative real). The maximum absolute error when computing the solution.
%                               If this variable is NaN, it has not been calculated (solution is not unique).
%
%                   BMatrix (NumEq x (NumEq*NumLag)) matrix with the B(i) matrices of reduced form coefficients.
%                            These will be used to set of the matrices for the state equation of the state-space
%                            representation.
%
%                   SMatrix (NumEq x (NumEq*(NumLag+1))) matrix with the S(i) matrices for the observed structure.
%                            S(0) will be used to create B(0), where B(0)*B(0)' = Q in the state space representation.
%
% OPTIONAL OUTPUT ModelEigenvalues (structure) with fields "Roots", "NumExact", "NumNumeric", "LargeRoots",
%                            "DimCompanion", "DimStability", "Amplitude", and "Period".
%
%
% ACKNOWLEDGEMENTS: This function is based on Volker Wieland's Matlab script "vwsolve.m"
%
%
%                       Written by: Anders Warne
%                                   New Area Wide Model Project
%                                   DG-R/EMO
%                                   European Central Bank (ECB)
%                                   Email: anders.warne@ecb.europa.eu
%                                   Copyright © 2006-2009 European Central Bank.
%
%                       First version: July 18, 2006.
%                        This version: February 13, 2009.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% LICENSE INFORMATION:
%
%      YADA is free software: you can redistribute it and/or modify
%      it under the terms of the GNU General Public License as published by
%      the Free Software Foundation, either version 3 of the License, or
%      (at your option) any later version.
%
%      This program is distributed in the hope that it will be useful,
%      but WITHOUT ANY WARRANTY; without even the implied warranty of
%      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%      GNU General Public License for more details.
%
%      You should have received a copy of the GNU General Public License
%      along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%      YADA is released under the GNU General Public License, Version 3,
%      29 June 2007 <http://www.gnu.org/licenses/>. The current release of
%      the program was last modified by the ECB on the "This version" date
%      above.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% CHANGELOG:
%
% * 19-07-2006: Updated function documentation.
%
% * 25-08-2006: Changed the function names for the AiM files.
%
% * 15-09-2006: Changed the upper bound for the eigenvalues from "1 + 1.e-6" to "1".
%
% * 19-09-2006: Added the optional output "ModelEigenvalues" with DSGE model properties.
%
% * 20-09-2006: Changed the upper bound for the eignevalues back to "1 + 1.e-6". That way, we
%               can allow for unit roots in the AiM solution, but remove them in the state
%               equation by not selecting an equation with a unit root.
%
% * 17-11-2006: Updated the function documentation.
%
% * 30-10-2007: Updated the documentation.
%
% * 16-05-2008: Added the "AIMTolerance" input variable.
%
% * 23-05-2008: Updated the documentation.
%
% * 27-05-2008: Added the mcode value 7 when a has infinite or NaN entries. This avoids
%               running the eig function on such a matrix.
%
% * 03-07-2008: All output variables are now initialized.
%
% * 05-12-2008: Updated the documentation.
%
% * 30-01-2009: Updated the documentation.
%
% * 13-02-2009: Made sure that the output matrices from running the "compute_aim_matrices"
%               function are always real. If not, mcode returns the values 81 and 82
%               respectively.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

mcode = 0;
MaxAbsError = NaN;
BMatrix = [];
SMatrix = [];
ModelEigenvalues = [];
%
% Numerical tolerances for AiM
%
% condn is the tolerance number used to compute the numerical part of the left invariant subspace
% uprbnd is an upper bound for the modulus of the roots of the reduced form solution
%
%condn  = 1.e-8;
condn = 10^(-AIMTolerance-2);
%uprbnd = 1 + 1.e-6;
uprbnd = 1+(10^(-AIMTolerance));
%
% Construct structural coefficient matrices.
%
[GMatrix,HMatrix] = compute_aim_matrices(ModelParameters);
if isreal(GMatrix)==0;
   mcode = 8;
   return;
end;
if isreal(HMatrix)==0;
   mcode = 8;
   return;
end;
%
[rh,ch] = size(HMatrix);
[rg,cg] = size(GMatrix);
cof = zeros(rh,ch);
cof(1:rg,1:cg) = GMatrix;
cof = cof + HMatrix;
%
% create coefficients (Bmatrix), roots (rts), dimension of companion matrix (ia),
% number of exact shiftrights (nex), number of numeric shiftrights (nnum),
% number of large roots (lgrts), and check uniqueness of solution (mcode)
%
% mcode = 1:      Unique solution
% mocde = 2:      Roots not correctly computed by real_schur
% mcode = 3:      Too many big roots
% mcode = 35:     Too many big roots, and q(:,right) is singular
% mcode = 4:      Too few big roots
% mcode = 45:     Too few big roots, and q(:,right) is singular
% mcode = 5:      q(:,right) is singular
% mcode = 61:     Too many exact shiftrights
% mcode = 62:     Too many numeric shiftrights
% mcode = 7:      The a matrix has infinite or NaN values. eigenvalues
%                 cannot be calculated.
% mcode = 8:      compute_aim_matrices returns complex values.
% else            Return code not properly specified
%
[BMatrix,rts,ia,nex,nnum,lgrts,mcode] = AiMEigenvalues(cof,NumEq,NumLag,NumLead,condn,uprbnd);
%
% check the accuracy of the solution
%
if mcode==1;
   MaxAbsError = CheckAiMSolution(NumEq,NumLag,NumLead,cof,BMatrix);
   %
   % Compute observable structure
   %
   SMatrix = AiMObservedStructure(cof,BMatrix,NumEq,NumLag,NumLead);
else;
   SMatrix = [];
end;

if nargout>4;
   ModelEigenvalues.Roots = rts;
   ModelEigenvalues.NumExact = nex;
   ModelEigenvalues.NumNumeric = nnum;
   ModelEigenvalues.LargeRoots = lgrts;
   ModelEigenvalues.DimCompanion = ia;
   ModelEigenvalues.DimStability = nex+nnum+lgrts-(NumEq*NumLead);
   %
   % calculate the amplitude and the periodicity
   %
%   epsi  = 1.e-8;
   epsi = 10^(-AIMTolerance-2);
   %
   lambda = rts;
   amp = abs(lambda);
   keep = find(lambda);
   lambda = lambda(keep);
   amp = amp(keep);
   %
   n = length(keep);
   angel = zeros(n,1);
   period = zeros(n,1)+NaN;
   %
   i = find(abs(lambda)>epsi);
   %
   for j = i
      angel(j) = abs(angle(lambda(j)));
   end
   %
   k = find(angel>epsi);
   if ~isempty(k);
      period(k) = (2*pi)./angel(k);
   end;
   ModelEigenvalues.Amplitude = amp;
   ModelEigenvalues.Period = period;
end;

%
% end of AiMSolver.m
%
