function [ys,check,penlt] = mze_steadystate(ys,exe)
% stephane [DOT] adjemian [AT] ens [DOT] fr
global M_ AutoregressiveProcesses

persistent idx idx_1 idx_2 idx_3 idx_4 idx_5 idx_6 idx_7 idx_8 idx_9 
persistent idx_10 idx_11 idx_12
persistent ida ida_table_in_1 ida_table_in_2 ida_table_out_1 ida_table_out_2 
persistent NumberOfShocks NumberOfParameters NumberOfEndogenousVariables
persistent FinalGoodLagrangeMultiplier_ss AveragedRelativePrices_ss PriceDistorsion_ss
persistent NablaPrice_ss OptimalRelativePrice_ss EmploymentAgencyLagrangeMultiplier_ss
persistent RealWageGrowthFactor_ss AveragedRelativeWages_ss OptimalRelativeRealWage_ss 
persistent WageDistorsion_ss NablaWage_ss TobinQ_ss OutputGap_ss InvestmentCost_ss 
persistent dInvestmentCost_ss EfficientInvestmentCost_ss dEfficientInvestmentCost_ss
persistent load_parameters fill_ys


if isempty(idx)
    idx_1 = strmatch('first_derivate_depreciation_rate_ss',M_.param_names,'exact');
    idx_2 = strmatch('beta',M_.param_names,'exact');
    idx_3 = strmatch('labour_supply_ss',M_.param_names,'exact');
    idx_4 = strmatch('ssGDP',M_.param_names,'exact');
    idx_5 = strmatch('thetaf',M_.param_names,'exact');
    idx_6 = strmatch('thetas',M_.param_names,'exact');
    idx_7 = strmatch('alpha',M_.param_names,'exact');
    idx_8 = strmatch('psif',M_.param_names,'exact');
    idx_9 = strmatch('psis',M_.param_names,'exact');
    idx_10 = strmatch('xip',M_.param_names,'exact');
    idx_11 = strmatch('xiw',M_.param_names,'exact');
    idx_12 = strmatch('steady_state_nominal_interest_factor',M_.param_names,'exact');
    if isempty(idx_4)
        idx_4 = 0;
    end
    NumberOfShocks = size(AutoregressiveProcesses,1);
    NumberOfParameters = M_.param_nbr;
    NumberOfEndogenousVariables = M_.endo_nbr;
    load_parameters = [];
    for i = 1:NumberOfParameters
        load_parameters = [ load_parameters deblank(M_.param_names(i,:)) ' = M_.params(' int2str(i) '); '];    
    end
    fill_ys = [];
    for i = 1:NumberOfEndogenousVariables
        fill_ys = [fill_ys 'ys(' int2str(i) ') = ' deblank(M_.endo_names(i,:)) '_ss' '; '];
    end
    SumOfDepreciationRates = 0 ;
    AveragedDepreciationRate = 0;
    Counter = 0;
end


% Do not call matlab's routine!
beta = NaN;
alpha = NaN;

eval(load_parameters);

ys = zeros(NumberOfEndogenousVariables,1);
check = 0;

% Compute the autoregressive parameters from the autocorrelation functions:
if isempty(ida)
    ida_table_in_1  = cell(NumberOfShocks,1);
    ida_table_in_2  = cell(NumberOfShocks,1);
    ida_table_out_1 = cell(NumberOfShocks,1);
    ida_table_out_2 = cell(NumberOfShocks,1);
end
for shock = 1:NumberOfShocks
    NumberOfLags = AutoregressiveProcesses{shock,3}; 
    Rho = NaN(NumberOfLags,1);
    for lag = 1:NumberOfLags
        if isempty(ida)
            tmp = strmatch([AutoregressiveProcesses{shock,2} '_rho' int2str(lag)],M_.param_names,'exact');
            ida_table_in_1(shock,lag) = {tmp};
        end
        Rho(lag) = M_.params(ida_table_in_1{shock,lag});
    end
    if isempty(ida)
        tmp = strmatch([AutoregressiveProcesses{shock,2} '_std'],M_.param_names,'exact');
        ida_table_in_2(shock,1) = {tmp};
    end
    Variance = M_.params(ida_table_in_2{shock,1})*M_.params(ida_table_in_2{shock,1});
    [ InnovationVariance, AutoregressiveParameters ] = autoregressive_process_specification(Variance,Rho,NumberOfLags);
    if InnovationVariance<0
        check = abs(InnovationVariance);
        % disp('SteadyState problem (specification of an autoregressive process)!')
        % InnovationVariance
        return
    end
    if isempty(ida)
        tmp = strmatch([AutoregressiveProcesses{shock,2} '_i_std'],M_.param_names,'exact');
        ida_table_out_1(shock,1) = {tmp};
    end
    M_.params(ida_table_out_1{shock,1}) = sqrt(InnovationVariance);
    for lag = 1:NumberOfLags
        if isempty(ida)
            tmp = strmatch([AutoregressiveProcesses{shock,2} '_ar' int2str(lag)],M_.param_names,'exact');
            ida_table_out_2(shock,lag) = {tmp};
        end
        M_.params(ida_table_out_2{shock,lag}) = AutoregressiveParameters(lag); 
    end
end
ida = 1;

steady_state_nominal_interest_factor = steady_state_real_interest_factor * inflation_target_ss ; 

thetaf = steady_state_mark_up_f/(steady_state_mark_up_f-1) ;
thetas = steady_state_mark_up_s/(steady_state_mark_up_s-1) ;

alpha = 1 - labour_tax_ss*steady_state_mark_up_f*steady_state_labour_share ;
if alpha>1-eps
    check = (alpha-1);
    disp('SteadyState problem (alpha is greater or equal to one)!')
    [alpha,labour_tax_ss,steady_state_mark_up_f, steady_state_labour_share]
    return
elseif alpha<eps
    check = abs(alpha);
    disp('SteadyState problem (alpha is negative or zero)!')
    [alpha,labour_tax_ss,steady_state_mark_up_f, steady_state_labour_share]
    return
end

psif = - kimball_curvature_f / thetaf ;
psis = - kimball_curvature_s / thetas ;

xip  = (price_contract_length - 1)/price_contract_length ;
xiw  = (wage_contract_length - 1)/wage_contract_length ;

if isempty(idx)
    FinalGoodLagrangeMultiplier_ss = 1 ;
    AveragedRelativePrices_ss = 1 ;
    PriceDistorsion_ss = 1 ;
    NablaPrice_ss = 1 ;
    OptimalRelativePrice_ss = 1 ;
    EmploymentAgencyLagrangeMultiplier_ss = 1 ;
    RealWageGrowthFactor_ss = 1;
    AveragedRelativeWages_ss = 1 ;
    OptimalRelativeRealWage_ss = 1 ;
    WageDistorsion_ss = 1 ;
    NablaWage_ss = 1 ;
    TobinQ_ss = 1;
    OutputGap_ss = 1 ;
    InvestmentCost_ss = 0;
    dInvestmentCost_ss = 0;
    EfficientInvestmentCost_ss = 0;
    dEfficientInvestmentCost_ss = 0;
    idx = 1;
end

RealMarginalCost_ss = 1/steady_state_mark_up_f ;
NonOptimizingFirmsPriceGrowth_ss = inflation_target_ss ;
InflationFactor_ss = inflation_target_ss ;
NonOptimizingUnionsWageGrowth_ss = inflation_target_ss ;

ssratio_HouseholdRealWage_RealGrossWage = 1/steady_state_mark_up_s ; 

NominalInterestFactor_ss = steady_state_nominal_interest_factor ;
beta = inflation_target_ss/(risk_premium_ss*NominalInterestFactor_ss)*(1+production_efficiency_growth_ss)^(sigmac);
if beta/(1+production_efficiency_growth_ss)^(sigmac)>1
    check = (beta/(1+production_efficiency_growth_ss)^(sigmac)-1);
    disp('SteadyState problem (the discount factor is greater than (1+g)^sigmac)!')
    [beta inflation_target_ss risk_premium_ss NominalInterestFactor_ss production_efficiency_growth_ss sigmac]
    return
end
    
RealInterestFactor_ss = NominalInterestFactor_ss / inflation_target_ss ;  

CapacityUtilizationFactor_ss = capacity_utilization_factor_ss ;

DepreciationRate_ss = depreciation_rate_ss;

CapitalReturnRate_ss = ((1+production_efficiency_growth_ss)^sigmac-beta*(1-depreciation_rate_ss)) ...
    /(beta*income_tax_ss*capacity_utilization_factor_ss);

dDepreciationRate_ss = CapitalReturnRate_ss*income_tax_ss;

RealGrossWage_ss = (RealMarginalCost_ss*(alpha/CapitalReturnRate_ss)^alpha)^(1/(1-alpha))*(1-alpha)/labour_tax_ss*production_efficiency_ss;

HouseholdRealWage_ss = ssratio_HouseholdRealWage_RealGrossWage * RealGrossWage_ss;

ssratio_CapitalDemand_HouseholdLabourSupply = labour_tax_ss*RealGrossWage_ss/CapitalReturnRate_ss*alpha/(1-alpha) ;

ssratio_GDP_HouseholdLabourSupply = ssratio_CapitalDemand_HouseholdLabourSupply^alpha * production_efficiency_ss^(1-alpha) ;

ssratio_Investment_GDP = (production_efficiency_growth_ss + depreciation_rate_ss) ...
    /capacity_utilization_factor_ss*ssratio_CapitalDemand_HouseholdLabourSupply^(1-alpha)/production_efficiency_ss^(1-alpha) ; 

if (ssratio_Investment_GDP > 1-public_spending_ss)
    check = ssratio_Investment_GDP - (1-public_spending_ss);
    disp('SteadyState problem (the investment ratio is greater or equal to one)!')
    [ssratio_Investment_GDP, public_spending_ss, production_efficiency_growth_ss, depreciation_rate_ss, capacity_utilization_factor_ss]
    return
elseif (ssratio_Investment_GDP < 0)
    check = -ssratio_Investment_GDP  ;
    disp('SteadyState problem (the investment ratio is negative or zero)!')
    [ssratio_Investment_GDP]
    return
end

ssratio_Consumption_GDP = 1 - public_spending_ss - ssratio_Investment_GDP;

labour_supply_ss = HouseholdRealWage_ss/(ssratio_Consumption_GDP*ssratio_GDP_HouseholdLabourSupply) ...
                         * income_tax_ss * labour_income_tax_ss / consumption_tax_ss ...
                         /(1-eta/(1+production_efficiency_growth_ss));

HouseholdLabourSupply_ss = household_labour_supply_ss;
LabourDemand_ss = HouseholdLabourSupply_ss;
EmploymentAgencyLabourSupply_ss = HouseholdLabourSupply_ss;

GDP_ss = ssratio_GDP_HouseholdLabourSupply * HouseholdLabourSupply_ss ;
Consumption_ss = ssratio_Consumption_GDP * GDP_ss ;
Investment_ss = ssratio_Investment_GDP * GDP_ss ;
HouseholdLagrangeMultiplier_ss = (1/consumption_tax_ss)*(Consumption_ss*(1-eta/(1+production_efficiency_growth_ss)))^(-sigmac) ...
    * exp(labour_supply_ss*(sigmac-1)/(1+sigmal));
CapitalDemand_ss = ssratio_CapitalDemand_HouseholdLabourSupply * HouseholdLabourSupply_ss ;
CapitalStock_ss = (1+production_efficiency_growth_ss)/(production_efficiency_growth_ss+depreciation_rate_ss)*Investment_ss ;

Z1_ss = RealMarginalCost_ss*HouseholdLagrangeMultiplier_ss*GDP_ss ...
        /(1-beta*xip*(1+production_efficiency_growth_ss)^(1-sigmac)) ;
Z2_ss = Z1_ss/RealMarginalCost_ss ;
Z3_ss = Z2_ss ;

H1_ss = HouseholdRealWage_ss*HouseholdLagrangeMultiplier_ss*EmploymentAgencyLabourSupply_ss ...
        /(1-beta*xiw*(1+production_efficiency_growth_ss)^(1-sigmac)) ;
H2_ss = H1_ss/ssratio_HouseholdRealWage_RealGrossWage ;
H3_ss = H2_ss ;

ProductionEfficiencyGrowth_ss = (1+production_efficiency_growth_ss) ;
ProductionEfficiencyCycle_ss = production_efficiency_ss ;
ConsumptionTax_ss = consumption_tax_ss ;
LabourIncomeTax_ss = labour_income_tax_ss ;
IncomeTax_ss = income_tax_ss ;
RiskPremium_ss = risk_premium_ss ;
LabourSupplyShock_ss = labour_supply_ss ;
InvestmentRelativePrice_ss = investment_relative_price_ss ;
InvestmentEfficiencyShock_ss = investment_efficiency_ss ;

LabourTax_ss = labour_tax_ss ;
PriceMarkUpShock_ss = price_mark_up_ss ;
WageMarkUpShock_ss = wage_mark_up_ss ;
InflationTarget_ss = inflation_target_ss ;
TaylorShock_ss = taylor_ss ;
PublicSpendingShare_ss = public_spending_ss ;
PublicSpending_ss = GDP_ss * PublicSpendingShare_ss ;

EfficientConsumption_ss = Consumption_ss ;
EfficientHouseholdLabourSupply_ss = HouseholdLabourSupply_ss ;
EfficientHouseholdLagrangeMultiplier_ss = HouseholdLagrangeMultiplier_ss ;
EfficientInvestment_ss = Investment_ss;
EfficientCapitalStock_ss = CapitalStock_ss;
EfficientDepreciationRate_ss = DepreciationRate_ss;
dEfficientDepreciationRate_ss = dDepreciationRate_ss;
EfficientCapacityUtilizationFactor_ss = CapacityUtilizationFactor_ss;
EfficientCapitalReturnRate_ss = CapitalReturnRate_ss;
EfficientHouseholdRealWage_ss = HouseholdRealWage_ss;
EfficientTobinQ_ss = TobinQ_ss;
EfficientRealInterestFactor_ss = RealInterestFactor_ss;
EfficientLabourDemand_ss = LabourDemand_ss;
EfficientRealMarginalCost_ss = RealMarginalCost_ss;
EfficientGDP_ss = GDP_ss;
EfficientCapitalDemand_ss = CapitalDemand_ss;
EfficientRealGrossWage_ss = RealGrossWage_ss; 

HabitShock_ss = 1.0000;
PreferenceShock_ss = 1.0000;
DynamicWelfare_ss = ((Consumption_ss*(1-eta*habit_ss))^(1-sigmac))/(1-sigmac)*exp(labour_supply_ss*(sigmac-1)/(1+sigmal))/(1-beta);
CounterFactualConsumption_ss = ((DynamicWelfare_ss*(1-beta)*(1-sigmac)/exp(labour_supply_ss*(sigmac-1)/(1+sigmal)))^(1/(1-sigmac)))/(1-eta*habit_ss); 


% Steady state level of the observed variables.
dY_ss = 100*log(1+production_efficiency_growth_ss);
dC_ss = 100*log(1+production_efficiency_growth_ss);
dI_ss = 100*log(1+production_efficiency_growth_ss);
dW_ss = 100*log(1+production_efficiency_growth_ss); 
LoggedPI_ss = 100*log(inflation_target_ss);
LoggedR_ss  = 100*log(NominalInterestFactor_ss);
TUC_ss  = 100*capacity_utilization_factor_ss;
Lobs_ss = 4*household_labour_supply_ss;

% Update some of the parameters.
M_.params(idx_1) = dDepreciationRate_ss;
M_.params(idx_2) = beta ;
M_.params(idx_3) = labour_supply_ss ;
if idx_4
    M_.params(idx_4) = GDP_ss ;
end
M_.params(idx_5) = thetaf ;
M_.params(idx_6) = thetas ;
M_.params(idx_7) = alpha ;
M_.params(idx_8) = psif ;
M_.params(idx_9) = psis ;
M_.params(idx_10) = xip ;
M_.params(idx_11) = xiw ;
M_.params(idx_12) = steady_state_nominal_interest_factor ;

% Fill vector ys (steady state values).
eval(fill_ys);
