Fixed steady state in rbcii.mod (steady state levels in the unconstrained regime were wrong).

tests/ep/rbcii.mod

@@ -60,38 +60,14 @@ model;
 
 end;
 
-steady_state_model;
-efficiency = 0;
-Efficiency = effstar*exp(efficiency-mean_preserving_spread(rho,sigma2));
-// Compute steady state ratios.
-Output_per_unit_of_Capital=((1/beta-1+delta)/alpha)^(1/(1-psi));
-Consumption_per_unit_of_Capital=Output_per_unit_of_Capital-delta;
-Labour_per_unit_of_Capital=(((Output_per_unit_of_Capital/Efficiency)^psi-alpha)/(1-alpha))^(1/psi);
-Output_per_unit_of_Labour=Output_per_unit_of_Capital/Labour_per_unit_of_Capital;
-Consumption_per_unit_of_Labour=Consumption_per_unit_of_Capital/Labour_per_unit_of_Capital;
-
-// Compute steady state share of capital.
-ShareOfCapital=alpha/(alpha+(1-alpha)*Labour_per_unit_of_Capital^psi);
-
-/// Compute steady state of the endogenous variables.
-Labour=1/(1+Consumption_per_unit_of_Labour/((1-alpha)*theta/(1-theta)*Output_per_unit_of_Labour^(1-psi)));
-Consumption = Consumption_per_unit_of_Labour*Labour;
-Capital = Labour/Labour_per_unit_of_Capital;
-Output = Output_per_unit_of_Capital*Capital;
-Investment = delta*Capital;
-ExpectedTerm = beta*((((Consumption^theta)*((1-Labour)^(1-theta)))^(1-tau))/Consumption)*(alpha*((Output/Capital)^(1-psi))+1-delta);
-Output1 = Output;
-Output2 = Output;
-Labour1 = Labour;
-Labour2 = Labour;
-Consumption1 = Consumption;
-Consumption2 = Consumption;
-end;
+// Write analytical steady state file (without globals)
+options_.steadystate_flag = 2;
+copyfile('rbcii_steady_state.m','rbcii_steadystate2.m');
 
 @#if extended_path_version
 
     shocks;
-    var EfficiencyInnovation = sigma2;
+    var EfficiencyInnovation = 1;
     end;
 
     steady(nocheck);

tests/ep/rbcii_steady_state.m

@@ -0,0 +1,88 @@
+function [ys_, params, info] = rbcii_steadystate2(ys_, exo_, params)
+     
+    % Flag initialization (equal to zero if the deterministic steady state exists) 
+    info = 0;
+    
+    % efficiency
+    ys_(13)=0;
+    
+    % Efficiency
+    ys_(12)=params(8);
+    
+    % Steady state ratios 
+    Output_per_unit_of_Capital=((1/params(1)-1+params(6))/params(4))^(1/(1-params(5)));
+    Consumption_per_unit_of_Capital=Output_per_unit_of_Capital-params(6);
+    Labour_per_unit_of_Capital=(((Output_per_unit_of_Capital/ys_(12))^params(5)-params(4))/(1-params(4)))^(1/params(5));
+    Output_per_unit_of_Labour=Output_per_unit_of_Capital/Labour_per_unit_of_Capital;
+    Consumption_per_unit_of_Labour=Consumption_per_unit_of_Capital/Labour_per_unit_of_Capital;
+
+    % Steady state share of capital revenues in total revenues (calibration check) 
+    ShareOfCapital=params(4)/(params(4)+(1-params(4))*Labour_per_unit_of_Capital^params(5));
+
+    % Steady state level of labour
+    ys_(3)=1/(1+Consumption_per_unit_of_Labour/((1-params(4))*params(2)/(1-params(2))*Output_per_unit_of_Labour^(1-params(5))));
+    
+    % Steady state level of consumption
+    ys_(4)=Consumption_per_unit_of_Labour*ys_(3);
+    
+    % Steady state level of physical capital stock
+    ys_(1)=ys_(3)/Labour_per_unit_of_Capital;
+    
+    % Steady state level of output
+    ys_(2)=Output_per_unit_of_Capital*ys_(1);
+    
+    % Steady state level of investment
+    ys_(5)=params(6)*ys_(1);
+    
+    % Steady state level of the expected term appearing in the Euler equation
+    ys_(14)=params(1)*(ys_(4)^params(2)*(1-ys_(3))^(1-params(2)))^(1-params(3))/ys_(4)*(1+params(4)*(ys_(2)/ys_(1))^(1-params(5))-params(6));
+
+    % Steady state level of output in the unconstrained regime (positive investment)
+    ys_(6)=ys_(2);
+
+    % Steady state level of labour in the unconstrained regime
+    ys_(7)=ys_(3);
+    
+    % Steady state level of consumption in the unconstrained regime 
+    ys_(8)=ys_(4);
+        
+    % Steady state level of labour in the constrained regime (noinvestment)
+    [lss,info] = l_solver(ys_(3),params(4),params(5),params(2),params(8),ys_(1),100);
+    if info, return, end
+    ys_(10) = lss;
+
+    % Steady state level of consumption in the constrained regime
+    ys_(11)=params(8)*(params(4)*ys_(1)^params(5)+(1-params(4))*ys_(10)^params(5))^(1/params(5));
+    
+    % Steady state level of output in the constrained regime
+    ys_(9)=ys_(11);
+
+end
+
+
+
+
+function r = p0(labour,alpha,psi,theta,effstar,kstar)
+    r = labour * ( alpha*kstar^psi/labour^psi + 1-alpha + theta*(1-alpha)/(1-theta)/effstar^psi ) - theta*(1-alpha)/(1-theta)/effstar^psi;
+end
+    
+function d = p1(labour,alpha,psi,theta,effstar,kstar)
+    d = alpha*(1-psi)*kstar^psi/labour^psi + 1-alpha + theta*(1-alpha)/(1-alpha)/effstar^psi;
+end
+
+function [labour,info] = l_solver(labour,alpha,psi,theta,effstar,kstar,maxiter)
+    iteration = 1; info = 0;
+    r = p0(labour,alpha,psi,theta,effstar,kstar);
+    condition = abs(r);
+    while condition
+        if iteration==maxiter
+            info = 1;
+            break
+        end
+        d = p1(labour,alpha,psi,theta,effstar,kstar);
+        labour = labour - r/d;
+        r = p0(labour,alpha,psi,theta,effstar,kstar);
+        condition = abs(r)>1e-9;
+        iteration = iteration + 1; 
+    end
+end