% write_ss_initval.m
% Run from BBCDL_model_excl2ndorder/.
% Computes the Ramsey steady state independently of Dynare's M_ structure
% and rewrites the initval block in the .mod file with true SS values.

addpath parameterfiles/

%% Parameters
paramfile_BB;
omega_welf1 = 0.5;
omega_welf2 = 1 - omega_welf1;
nbeta       = c1beta;

%% Non-lmult SS values
% Variables that appear as exp(var) in the .mod file → stored in log form.
% Variables that appear directly                     → stored as levels.

c2rerss     = 1 / c1rerss;

iv_Util1       = Util1ss;
iv_Util2       = Util2ss;
iv_Util        = omega_welf1*Util1ss + omega_welf2*Util2ss;
iv_Welf1       = Util1ss / (1-nbeta);
iv_Welf2       = Util2ss / (1-nbeta);
iv_Welf        = iv_Util   / (1-nbeta);
iv_Welfcon1    = Welfcon1;
iv_instr1      = c1pidss;       % level = 1
iv_instr2      = c2pidss;       % level = 1

iv_c1c         = log(c1css);
iv_c1cd        = log(c1cdss);
iv_c1mc        = log(c1mcss);
iv_c1y         = log(c1yss);
iv_c1g         = log(c1gss);
iv_c1rpcpd     = log(c1rpcpdss);   % log(1) = 0
iv_c1pid       = c1pidss;          % level = 1
iv_c1hp        = log(c1hpss);
iv_c1gp        = log(c1gpss);
iv_c1pstar     = log(c1pstarss);   % log(1) = 0
iv_c1deltad    = log(c1deltadss);  % log(1) = 0
iv_c1rn        = c1rnss;           % level
iv_c1rer       = log(c1rerss);     % log(1) = 0
iv_c1dcore     = c1pistar;         % level = 1
iv_c1ngdp      = log(c1yss);
iv_c1pcore     = 0;                % log(1) = 0
iv_c1dnxf      = 0;                % log(1) = 0

iv_c2c         = log(c2css);
iv_c2cd        = log(c2cdss);
iv_c2mc        = log(c2mcss);
iv_c2y         = log(c2yss);
iv_c2g         = log(c2gss);
iv_c2rpcpd     = log(c2rpcpdss);   % log(1) = 0
iv_c2pid       = c2pidss;          % level = 1
iv_c2hp        = log(c2hpss);
iv_c2gp        = log(c2gpss);
iv_c2pstar     = 0;                % log(1) = 0
iv_c2deltad    = log(c2deltadss);  % log(1) = 0
iv_c2rn        = c2rnss;           % level
iv_c2rer       = log(c2rerss);     % log(1/1) = 0
iv_c2dcore     = c2pistar;         % level = 1
iv_c2ngdp      = log(c2yss);
iv_c2pcore     = 0;
iv_c2dnxf      = 0;

% Lag variables equal current period values at SS
iv_c1deltadlag = iv_c1deltad;
iv_c1ngdplag   = iv_c1ngdp;
iv_c1pcorelag  = iv_c1pcore;
iv_c1rerlag    = iv_c1rer;
iv_c1rnlag     = iv_c1rn;
iv_c1rpcpdlag  = iv_c1rpcpd;
iv_c1ylag      = iv_c1y;
iv_c2deltadlag = iv_c2deltad;
iv_c2ngdplag   = iv_c2ngdp;
iv_c2pcorelag  = iv_c2pcore;
iv_c2rerlag    = iv_c2rer;
iv_c2rnlag     = iv_c2rn;
iv_c2rpcpdlag  = iv_c2rpcpd;
iv_c2ylag      = iv_c2y;

% Exogenous shocks = 0 at SS
iv_c1shockgy   = 0;
iv_c1z         = 0;
iv_c1taupvar   = 0;
iv_c2shockgy   = 0;
iv_c2z         = 0;
iv_c2taupvar   = 0;

%% Set _ss variables required by lmss.m
c1c_ss = iv_c1c;    c1cd_ss = iv_c1cd;    c1mc_ss = iv_c1mc;
c1y_ss = iv_c1y;    c1g_ss  = iv_c1g;     c1rpcpd_ss = iv_c1rpcpd;
c1hp_ss = iv_c1hp;  c1gp_ss = iv_c1gp;    c1pstar_ss = iv_c1pstar;
c1pid_ss = iv_c1pid; c1deltad_ss = iv_c1deltad; c1rn_ss = iv_c1rn;
c1rer_ss = iv_c1rer; c1dcore_ss = iv_c1dcore;
c1ngdp_ss = iv_c1ngdp; c1pcore_ss = iv_c1pcore; c1dnxf_ss = iv_c1dnxf;
c1z_ss = iv_c1z; c1shockgy_ss = iv_c1shockgy; c1taupvar_ss = iv_c1taupvar;

c2c_ss = iv_c2c;    c2cd_ss = iv_c2cd;    c2mc_ss = iv_c2mc;
c2y_ss = iv_c2y;    c2g_ss  = iv_c2g;     c2rpcpd_ss = iv_c2rpcpd;
c2hp_ss = iv_c2hp;  c2gp_ss = iv_c2gp;    c2pstar_ss = iv_c2pstar;
c2pid_ss = iv_c2pid; c2deltad_ss = iv_c2deltad; c2rn_ss = iv_c2rn;
c2rer_ss = iv_c2rer; c2dcore_ss = iv_c2dcore;
c2ngdp_ss = iv_c2ngdp; c2pcore_ss = iv_c2pcore; c2dnxf_ss = iv_c2dnxf;
c2z_ss = iv_c2z; c2shockgy_ss = iv_c2shockgy; c2taupvar_ss = iv_c2taupvar;

Util1_ss = iv_Util1; Util2_ss = iv_Util2;
Util_ss  = iv_Util;  Welf_ss  = iv_Welf;
Welf1_ss = iv_Welf1; Welf2_ss = iv_Welf2;
Welfcon1_ss = iv_Welfcon1;
instr1_ss = iv_instr1; instr2_ss = iv_instr2;
c1deltadlag_ss = iv_c1deltadlag; c1ngdplag_ss = iv_c1ngdplag;
c1pcorelag_ss  = iv_c1pcorelag;  c1rerlag_ss   = iv_c1rerlag;
c1rnlag_ss     = iv_c1rnlag;     c1rpcpdlag_ss = iv_c1rpcpdlag;
c1ylag_ss      = iv_c1ylag;
c2deltadlag_ss = iv_c2deltadlag; c2ngdplag_ss = iv_c2ngdplag;
c2pcorelag_ss  = iv_c2pcorelag;  c2rerlag_ss   = iv_c2rerlag;
c2rnlag_ss     = iv_c2rnlag;     c2rpcpdlag_ss = iv_c2rpcpdlag;
c2ylag_ss      = iv_c2ylag;

%% Solve for lmult SS values
BBCDLmodelcomp_ramsey_c1pid_c2pid_lmss;
% lmss.m ends with:  lmss_vals = -lmss_mat \ lmss_vec;
% and assigns lmult1 ... lmult51 from lmss_vals.

%% Build initval block string
fv = @(x) sprintf('%.15g', x);

lines = { ...
    'initval;', ...
    ['Util = '       fv(iv_Util)       ';'], ...
    ['Welf = '       fv(iv_Welf)       ';'], ...
    ['instr1 = '     fv(iv_instr1)     ';'], ...
    ['instr2 = '     fv(iv_instr2)     ';'], ...
    ['Util1 = '      fv(iv_Util1)      ';'], ...
    ['Util2 = '      fv(iv_Util2)      ';'], ...
    ['c1c = '        fv(iv_c1c)        ';'], ...
    ['c1cd = '       fv(iv_c1cd)       ';'], ...
    ['c1mc = '       fv(iv_c1mc)       ';'], ...
    ['c1y = '        fv(iv_c1y)        ';'], ...
    ['c1g = '        fv(iv_c1g)        ';'], ...
    ['c1rpcpd = '    fv(iv_c1rpcpd)    ';'], ...
    ['c1pid = '      fv(iv_c1pid)      ';'], ...
    ['c1hp = '       fv(iv_c1hp)       ';'], ...
    ['c1gp = '       fv(iv_c1gp)       ';'], ...
    ['c1pstar = '    fv(iv_c1pstar)    ';'], ...
    ['c1deltad = '   fv(iv_c1deltad)   ';'], ...
    ['c1rn = '       fv(iv_c1rn)       ';'], ...
    ['c1rer = '      fv(iv_c1rer)      ';'], ...
    ['c1dcore = '    fv(iv_c1dcore)    ';'], ...
    ['c1ngdp = '     fv(iv_c1ngdp)     ';'], ...
    ['c1pcore = '    fv(iv_c1pcore)    ';'], ...
    ['c1dnxf = '     fv(iv_c1dnxf)     ';'], ...
    ['c2c = '        fv(iv_c2c)        ';'], ...
    ['c2cd = '       fv(iv_c2cd)       ';'], ...
    ['c2mc = '       fv(iv_c2mc)       ';'], ...
    ['c2y = '        fv(iv_c2y)        ';'], ...
    ['c2g = '        fv(iv_c2g)        ';'], ...
    ['c2rpcpd = '    fv(iv_c2rpcpd)    ';'], ...
    ['c2pid = '      fv(iv_c2pid)      ';'], ...
    ['c2hp = '       fv(iv_c2hp)       ';'], ...
    ['c2gp = '       fv(iv_c2gp)       ';'], ...
    ['c2pstar = '    fv(iv_c2pstar)    ';'], ...
    ['c2deltad = '   fv(iv_c2deltad)   ';'], ...
    ['c2rn = '       fv(iv_c2rn)       ';'], ...
    ['c2rer = '      fv(iv_c2rer)      ';'], ...
    ['c2dcore = '    fv(iv_c2dcore)    ';'], ...
    ['c2ngdp = '     fv(iv_c2ngdp)     ';'], ...
    ['c2pcore = '    fv(iv_c2pcore)    ';'], ...
    ['c2dnxf = '     fv(iv_c2dnxf)     ';'], ...
    ['Welfcon1 = '   fv(iv_Welfcon1)   ';'], ...
    ['c1deltadlag = ' fv(iv_c1deltadlag) ';'], ...
    ['c1ngdplag = '  fv(iv_c1ngdplag)  ';'], ...
    ['c1pcorelag = ' fv(iv_c1pcorelag) ';'], ...
    ['c1rerlag = '   fv(iv_c1rerlag)   ';'], ...
    ['c1rnlag = '    fv(iv_c1rnlag)    ';'], ...
    ['c1rpcpdlag = ' fv(iv_c1rpcpdlag) ';'], ...
    ['c1ylag = '     fv(iv_c1ylag)     ';'], ...
    ['c2deltadlag = ' fv(iv_c2deltadlag) ';'], ...
    ['c2ngdplag = '  fv(iv_c2ngdplag)  ';'], ...
    ['c2pcorelag = ' fv(iv_c2pcorelag) ';'], ...
    ['c2rerlag = '   fv(iv_c2rerlag)   ';'], ...
    ['c2rnlag = '    fv(iv_c2rnlag)    ';'], ...
    ['c2rpcpdlag = ' fv(iv_c2rpcpdlag) ';'], ...
    ['c2ylag = '     fv(iv_c2ylag)     ';'], ...
    ['c1shockgy = '  fv(iv_c1shockgy)  ';'], ...
    ['c1z = '        fv(iv_c1z)        ';'], ...
    ['c1taupvar = '  fv(iv_c1taupvar)  ';'], ...
    ['c2shockgy = '  fv(iv_c2shockgy)  ';'], ...
    ['c2z = '        fv(iv_c2z)        ';'], ...
    ['c2taupvar = '  fv(iv_c2taupvar)  ';'], ...
};

lmult_lines = arrayfun(@(k) ['lmult' num2str(k) ' = ' fv(lmss_vals(k)) ';'], 1:51, 'UniformOutput', false);
lines = [lines, lmult_lines, {['Welf1 = ' fv(iv_Welf1) ';'], ['Welf2 = ' fv(iv_Welf2) ';'], 'end;'}];

new_block = strjoin(lines, newline);

%% Replace initval block in .mod file
modfile = 'BBCDLmodelcomp_ramsey_c1pid_c2pid.mod';
txt = fileread(modfile);

% Match from 'initval;' to the next 'end;'
txt_new = regexprep(txt, 'initval;.*?end;', new_block, 'dotall');

if strcmp(txt, txt_new)
    warning('initval block not found or unchanged — check the .mod file.');
else
    fid = fopen(modfile, 'w');
    fwrite(fid, uint8(txt_new));
    fclose(fid);
    fprintf('initval block updated in %s\n', modfile);
    fprintf('lmult range: [%.4g, %.4g]\n', min(lmss_vals), max(lmss_vals));
end
