% run.m:  Model: Sticky price model of the business cycle similar to the
%         setting developed in Ireland (2003).
%
%         Randomly-chosen starting values are used to find the
%         maximum of the model's likelihood function.
%
%         Output:
%         startingval
%         paramvec
%         likeval
%         errormat
%         sevecstarmat
%
%         Country: France
%
% O. Roehe (2010)


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Housekeeping
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%clear

%clc

%%
  global ct it mt pt rt bigt scalinv

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load Slovenia_centered2;

  %ct = Slovenia_centered2(:,1);
  %it = Slovenia_centered2(:,2);
  %mt = Slovenia_centered2(:,3);
 %pt = Slovenia_centered2(:,4);
%rt = Slovenia_centered2(:,5);

  bigt = length(ct);

  trend = 1:bigt;


  fstarmat = zeros(2500,1);     % Contains the log-likelihood values for
                                % each run

  tostarmat = zeros(20,2500); % Contains the transformed starting values
                                % for each run

  tstarmat = zeros(20,2500);    % Contains the paramerter estimates
                                % for each run

  flagmat = zeros(2500,1);     % Saves MATLAB error messages

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Generating starting Values (transformed parameters)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  for t0 = 1:2500

    t0

    theto = rand(10,1); % Generating starting values for gamma, phi_P,
                        % omega_mu, omega_pi, omega_y, sigma_a, sigma_e,
                        % Sigma_x, sigma_z, and sigma_v

    thetoscale = [ 1 60 1 3 0 0.001 0.001 0.1 0.1 0.001]';

    theto = thetoscale.*theto;

    betatr = 0.1094;

    gammatr = sqrt(theto(1)/(1-theto(1)));

    alphatr = 0.29;

    etatr = 1.5;

    thetatr = 5;

    gtr = 0;

    deltatr = sqrt(0.025/0.975);

    phiptr = theto(2)/100;

    phiktr = 0.35;

    mutr = 0.009;

    omegamtr = theto(3);

    omegaptr = theto(4);

    omegaytr = theto(5);

    etr = 4.4423;   % Set to match with re-meaned HP-filtered data (See DeJong and Dave, 2007)

    ztr = 0.4242;   % Set to match with re-meaned HP-filtered data (See DeJong and Dave, 2007)

    rhoatr = sqrt(0.9/0.1)/10;

    rhoetr = sqrt(0.9/0.1)/10;

    rhoxtr = sqrt(0.9/0.1)/10;

    rhoztr = sqrt(0.9/0.1)/10;

 	rhovtr = sqrt(0.9/0.1)/10;

    sigatr = theto(6);

    sigetr = theto(7);

    sigxtr = theto(8);

    sigztr = theto(9);

    sigvtr = theto(10);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Maximizes the (minimizes the negative)log likelihood function
% for the model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  bigtheto = [ betatr gammatr ...
               phiptr phiktr ...
               mutr omegamtr omegaptr omegaytr ...
               etr ztr ...
               rhoatr rhoetr rhoxtr rhoztr rhovtr ...
               sigatr sigetr sigxtr sigztr sigvtr]';

    startingval(:,t0) = bigtheto;

    save startingval startingval

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Chose between optimization algorithms to maximize
% the objective function:
%  a) fminunc (MATLAB)
%  or
%  b) csminwel (Chris Sims)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% a) fminunc - find the minimum of an unconstrained multivariable function
% (derivative based)

 options = optimset('Display','iter','LargeScale','off','MaxFunEvals',10000,'MaxIter',10000);

 [thetstar,fstar,exitflag] = fminunc(@llfn,bigtheto,options);

% b) Hybrid optimization algorithm developed by Chris Sims
% combining derivative based BFGS-method with a simplex algorithm :

%   crit = 1e-7;
%   nit = 150; %maximal number of iteration in CSMINWEL
%   H0=1e-4*eye(20);

%   [fstar,thetstar,gh,H,itct,fcount,retcodeh]=csminwel('llfn',bigtheto,H0,[],crit,nit);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Untransform the estimated parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  thetstar = real(thetstar);

  beta = (100*thetstar(1))^2/(1+(100*thetstar(1))^2);
  gamma =thetstar(2)^2/(1+thetstar(2)^2);
  phip = abs(thetstar(3))*100;
  phik = abs(thetstar(4))*100;
  mu = 1 + abs(thetstar(5));
  omegam = thetstar(6);
  omegap = thetstar(7);
  omegay = thetstar(8);
  e = abs(thetstar(9));
  z = 10000*abs(thetstar(10));
  rhoa = (10*thetstar(11))^2/(1+(10*thetstar(11))^2);
  rhoe = (10*thetstar(12))^2/(1+(10*thetstar(12))^2);
  rhox = (10*thetstar(13))^2/(1+(10*thetstar(13))^2);
  rhoz = (10*thetstar(14))^2/(1+(10*thetstar(14))^2);
  rhov = (10*thetstar(15))^2/(1+(10*thetstar(15))^2);
  siga = abs(thetstar(16));
  sige = abs(thetstar(17));
  sigx = abs(thetstar(18));
  sigz = abs(thetstar(19));
  sigv = abs(thetstar(20));

  tstar = [ beta gamma ...
            phip phik ...
            mu omegam omegap omegay ...
            e z ...
            rhoa rhoe rhox rhoz rhov ...
            siga sige sigx sigz sigv]';

  scalvec = [ 1 10 ...
              0.01 0.01 ...
              1 0.1 0.1 0.1 ...
              0.1 0.0001 ...
              1 1 1 1 1 ...
              10 10 10 10 10]';

  scalinv = inv(diag(scalvec));

  tstars = diag(scalvec)*tstar;

  fstar = llfnses(tstars);

  eee = 1e-5;

  epsmat = eee*eye(20);

  hessvec = zeros(20,1);

  for i = 1:20

    hessvec(i) = llfnses(tstars+epsmat(:,i));

  end

  hessmatl = zeros(20,20);

  for i = 1:20

    for j = 1:i

      hessmatl(i,j) = (llfnses(tstars+epsmat(:,i)+epsmat(:,j)) ...
                        -hessvec(i)-hessvec(j)+fstar)/(eee^2);

    end

  end

  hessmatu = hessmatl' - diag(diag(hessmatl));

  hessmatf = hessmatl + hessmatu;

  bighx = scalinv*inv(hessmatf)*scalinv';

  sevec = sqrt(diag(bighx));


   tstarmat(:,t0) = tstar;

   fstarmat(t0) = fstar;

   flagmat(t0) = exitflag;

   sevecstarmat(:,t0) = sevec;

   save tstarmat tstarmat
   save fstarmat fstarmat
   save flagmat flagmat
   save sevecstarmat sevecstarmat

  end