OLGModel14: tighten GE constraints

OLGModel14.m

@@ -28,27 +28,31 @@
 % Lets model agents from age 20 to age 100, so 81 periods
 Params.agejshifter=19; % Age 20 minus one. Makes keeping track of actual age easy in terms of model age
 Params.J=100-Params.agejshifter; % =81, Number of period in life-cycle
-n_d.household=51; % Endogenous labour choice (fraction of time worked)
-n_a.household=501; % Endogenous share holdings
+n_d.household=101; % Endogenous labour choice (fraction of time worked)
+n_a.household=201; % Endogenous share holdings
 % Exogenous labor productivity units shocks (next two lines)
 n_z.household=15; % AR(1) with age-dependent params
 vfoptions.n_e.household=3; % iid
 N_j.household=Params.J; % Number of periods in finite horizon
 
 % Grids to use for firm
 n_d.firm=101; % Dividend payment
-n_a.firm=501; % Capital holdings
+n_a.firm=201; % Capital holdings
 n_z.firm=11; % Productivity shock
 N_j.firm=Inf; % Infinite horizon
 
+%% Global parameters (applies to household and firm)
+% Risk-free rate of return
+Params.r=0.05;
+
 %% Parameters for household
 
 % Discount rate
-Params.beta = 1.01; % Changed to get S to increase nearer to 1 given r=0.05 (ran it with beta=0.99, got S=0.3, so increased this; note that it interacts with sj to give the actual discount factor)
+Params.beta = 0.95; % Changed to get S to increase nearer to 1 given r=0.05 (ran it with beta=0.99, got S=0.3, so increased this; note that it interacts with sj to give the actual discount factor)
 % Preferences
 Params.sigma = 2; % Coeff of relative risk aversion (curvature of consumption)
 Params.eta = 1.5; % Curvature of leisure (This will end up being 1/Frisch elasty)
-Params.psi = 10; % Weight on leisure
+Params.psi = 2; % Weight on leisure
 
 % Demographics
 Params.agej=1:1:Params.J; % Is a vector of all the agej: 1,2,3,...,J
@@ -106,7 +110,7 @@
 Params.capadjconstant=1.21; % term in the capital adjustment cost
 % Tax
 Params.tau_corp=0.34; % Tax rate on corporate earnings
-Params.phi=0.52; % Fraction of capital adjustment costs that can be deducted from corporate earnings
+Params.phi=0.5; % Fraction of capital adjustment costs that can be deducted from corporate earnings
 Params.tau_d=0.2; % Tax rate on dividends
 Params.tau_cg=0.2; % Tax rate on capital gains
 % Idiosyncatic productivity shocks
@@ -124,22 +128,24 @@
 
 % Some initial values/guesses for variables that will be determined in general eqm
 Params.pension=0.4; % Initial guess (this will be determined in general eqm)
-Params.r=0.05;
-Params.w=1;
-Params.AccidentBeq=0.03; % Accidental bequests (this is the lump sum transfer)
+Params.w=1; % Wages, determines (household) labor supply and (firm) demand
+Params.AccidentBeq=0.02; % Accidental bequests (this is the lump sum transfer)
 Params.G=0.1; % Government expenditure
 Params.firmbeta=1/(1+Params.r/(1-Params.tau_cg)); % 1/(1+r) but returns net of capital gains tax
-Params.D=0.2; % Dividends
+Params.D=0.2; % Dividends rate expected/received by households
 Params.P0=1;
-Params.Lhscale=0.3; % Scaling the household labor supply
+Params.Lhscale=0.21; % Scaling the household labor supply
 
 %% Grids for household
 
 % Grid for labour choice
 d_grid.household=linspace(0,1,n_d.household)'; % Notice that it is imposing the 0<=h<=1 condition implicitly
 
 % Grid for share holdings
-a_grid.household=10*(linspace(0,1,n_a.household).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
+a_grid_cubed=linspace(0,1,ceil(n_a.household/2)).^3; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
+a_grid_linear=linspace(1,10,floor(n_a.household/2)+1);
+a_grid.household=[a_grid_cubed, a_grid_linear(2:end)]';
+% a_grid.household=linspace(0,6,n_a.household)';
 
 % First, z, the AR(1) with age-dependent parameters
 [z_grid_J, pi_z_J] = discretizeLifeCycleAR1_FellaGallipoliPan(Params.rho_z,Params.sigma_epsilon_z,n_z.household,Params.J);
@@ -165,7 +171,9 @@
 
 %% Grids for firm
 d_grid.firm=linspace(0,1,n_d.firm)'; % Notice that it is imposing the d>=0 condition implicitly
-a_grid.firm=10*(linspace(0,1,n_a.firm).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
+a_grid_cubed=linspace(0,1,ceil(n_a.firm/2)).^3; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
+a_grid_linear=linspace(1,10,floor(n_a.firm/2)+1);
+a_grid.firm=[a_grid_cubed, a_grid_linear(2:end)]';
 
 [z_grid.firm,pi_z.firm] = discretizeAR1_FarmerToda(0,Params.rho_z_firm,Params.sigma_z_e_firm,n_z.firm);
 z_grid.firm=exp(z_grid.firm);
@@ -176,14 +184,14 @@
 % For households
 DiscountFactorParamNames.household={'beta','sj'};
 % Notice we use 'OLGModel14_HouseholdReturnFn'
-ReturnFn.household=@(h,sprime,s,z,e,sigma,psi,eta,agej,Jr,J,pension,r,w,P0,D,kappa_j,warmglow1,warmglow2,AccidentBeq,tau_l,tau_d,tau_cg,Lhscale)...
-    OLGModel14_HouseholdReturnFn(h,sprime,s,z,e,sigma,psi,eta,agej,Jr,J,pension,r,w,P0,D,kappa_j,warmglow1,warmglow2,AccidentBeq,tau_l,tau_d,tau_cg,Lhscale);
+ReturnFn.household=@(h,sprime,s,z,e,sigma,psi,eta,agej,Jr,J,pension,w,P0,D,kappa_j,warmglow1,warmglow2,AccidentBeq,r,tau_l,tau_d,tau_cg)...
+    OLGModel14_HouseholdReturnFn(h,sprime,s,z,e,sigma,psi,eta,agej,Jr,J,pension,w,P0,D,kappa_j,warmglow1,warmglow2,AccidentBeq,r,tau_l,tau_d,tau_cg);
 
 % For firms
 DiscountFactorParamNames.firm={'firmbeta'};
 % Notice we use 'OLGModel14_FirmReturnFn'
-ReturnFn.firm=@(d,kprime,k,z,w,delta,alpha_k,alpha_l,capadjconstant,tau_corp,phi,tau_d,tau_cg)...
-    OLGModel14_FirmReturnFn(d,kprime,k,z,w,delta,alpha_k,alpha_l,capadjconstant,tau_corp,phi,tau_d,tau_cg);
+ReturnFn.firm=@(d,kprime,k,z,w,D,delta,alpha_k,alpha_l,capadjconstant,tau_corp,phi,tau_d,tau_cg)...
+    OLGModel14_FirmReturnFn(d,kprime,k,z,w,D,delta,alpha_k,alpha_l,capadjconstant,tau_corp,phi,tau_d,tau_cg);
 
 %% Now solve the value function iteration problem, just to check that things are working before we go to General Equilbrium
 disp('Test ValueFnIter')
@@ -216,7 +224,7 @@
 StationaryDist=StationaryDist_Case1_FHorz_PType(jequaloneDist,AgeWeightsParamNames,PTypeDistParamNames,Policy,n_d,n_a,n_z,N_j,Names_i,pi_z,Params,simoptions);
 
 %% General eqm variables
-GEPriceParamNames={'r','pension','AccidentBeq','G','w','firmbeta','D','P0','Lhscale'}; 
+GEPriceParamNames={'pension','AccidentBeq','G','w','firmbeta','D','P0'}; 
 % We don't need P
 % We can get P from the equation that defines r as the return to the mutual fund
 % 1+r = (P0 +(1-tau_d)D - tau_cg(P0-P))/Plag
@@ -280,15 +288,13 @@
 p_eqm=HeteroAgentStationaryEqm_Case1_FHorz_PType(n_d, n_a, n_z, N_j, Names_i, [], pi_z, d_grid, a_grid, z_grid,jequaloneDist, ReturnFn, FnsToEvaluate, GeneralEqmEqns, Params, DiscountFactorParamNames, AgeWeightsParamNames, PTypeDistParamNames, GEPriceParamNames,heteroagentoptions, simoptions, vfoptions);
 % p_eqm contains the general equilibrium parameter values
 % Put this into Params so we can calculate things about the initial equilibrium
-Params.r=p_eqm.r;
 Params.pension=p_eqm.pension;
 Params.AccidentBeq=p_eqm.AccidentBeq;
 Params.G=p_eqm.G;
 Params.w=p_eqm.w;
 Params.firmbeta=p_eqm.firmbeta;
 Params.D=p_eqm.D;
 Params.P0=p_eqm.P0;
-Params.Lhscale=p_eqm.Lhscale;
 
 % Calculate a few things related to the general equilibrium.
 [V, Policy]=ValueFnIter_Case1_FHorz_PType(n_d,n_a,n_z,N_j, Names_i, d_grid, a_grid, z_grid, pi_z, ReturnFn, Params, DiscountFactorParamNames, vfoptions);
@@ -308,8 +314,8 @@
 %% Calculate some aggregates and print findings about them
 
 % Add consumption to the FnsToEvaluate
-FnsToEvaluate.Consumption.household=@(h,sprime,s,z,e,agej,Jr,pension,r,w,P0,D,kappa_j,AccidentBeq,tau_l,tau_d,tau_cg,Lhscale) OLGModel14_HouseholdConsumptionFn(h,sprime,s,z,e,agej,Jr,pension,r,w,P0,D,kappa_j,AccidentBeq,tau_l,tau_d,tau_cg,Lhscale);
-FnsToEvaluate.Income.household=@(h,sprime,s,z,e,agej,Jr,pension,r,w,P0,D,kappa_j,AccidentBeq,tau_l,tau_d,tau_cg,Lhscale) OLGModel14_HouseholdIncomeFn(h,sprime,s,z,e,agej,Jr,pension,r,w,P0,D,kappa_j,AccidentBeq,tau_l,tau_d,tau_cg,Lhscale);
+FnsToEvaluate.Consumption.household=@(h,sprime,s,z,e,agej,Jr,pension,w,P0,D,kappa_j,AccidentBeq,r,tau_l,tau_d,tau_cg) OLGModel14_HouseholdConsumptionFn(h,sprime,s,z,e,agej,Jr,pension,w,P0,D,kappa_j,AccidentBeq,r,tau_l,tau_d,tau_cg);
+FnsToEvaluate.Income.household=@(h,sprime,s,z,e,agej,Jr,pension,w,P0,D,kappa_j,AccidentBeq,r,tau_l,tau_d,tau_cg) OLGModel14_HouseholdIncomeFn(h,sprime,s,z,e,agej,Jr,pension,w,P0,D,kappa_j,AccidentBeq,r,tau_l,tau_d,tau_cg);
 
 AggVars=EvalFnOnAgentDist_AggVars_FHorz_Case1_PType(StationaryDist, Policy, FnsToEvaluate, Params, n_d, n_a, n_z,N_j, Names_i, d_grid, a_grid, z_grid,simoptions);
 

OLGModel14_FirmReturnFn.m

@@ -1,4 +1,4 @@
-function F=OLGModel14_FirmReturnFn(dividend,kprime,k,z,w,delta,alpha_k,alpha_l,capadjconstant,tau_corp,phi,tau_d,tau_cg)
+function F=OLGModel14_FirmReturnFn(dividend,kprime,k,z,w,D,delta,alpha_k,alpha_l,capadjconstant,tau_corp,phi,tau_d,tau_cg)
 % Whether we set it up so that dividends or equity issuance is the decision
 % variable is unimportant, here I use dividends as the decision variable.
 
@@ -32,6 +32,8 @@
 % Firms per-period objective
 if s>=0 % enforce that 'no share repurchases allowed'
     F=((1-tau_d)/(1-tau_cg))*dividend-s;
+    % Disfavor discrepencies between dividends paid and expected (D)
+    F=F-(D-dividend)^2;
 end
 
 % Note: dividend payments cannot be negative is enforced by the grid on

OLGModel14_HouseholdConsumptionFn.m

@@ -1,4 +1,4 @@
-function c=OLGModel14_HouseholdConsumptionFn(h,sprime,s,z,e,agej,Jr,pension,r,w,P0,D,kappa_j,AccidentBeq,tau_l,tau_d,tau_cg,Lhscale)
+function c=OLGModel14_HouseholdConsumptionFn(h,sprime,s,z,e,agej,Jr,pension,w,P0,D,kappa_j,AccidentBeq,r,tau_l,tau_d,tau_cg)
 % Replace assets with 'share holdings'
 % Get rid of progressive taxes
 % Add Lhnormalize
@@ -14,7 +14,7 @@
 
 if agej<Jr % If working age
     %consumption = labor income + accidental bequest + share holdings (including dividend) - capital gains tax - next period share holdings
-    c=(1-tau_l)*w*kappa_j*exp(z+e)*Lhscale*h+((1-tau_d)*D+P0)*(s+AccidentBeq) -tau_cg*(P0-Plag)*(s+AccidentBeq)-P*sprime; 
+    c=(1-tau_l)*w*kappa_j*exp(z+e)*h+((1-tau_d)*D+P0)*(s+AccidentBeq) -tau_cg*(P0-Plag)*(s+AccidentBeq)-P*sprime; 
 else % Retirement
     c=pension+((1-tau_d)*D+P0)*(s+AccidentBeq) -tau_cg*(P0-Plag)*(s+AccidentBeq) - P*sprime;
 end

OLGModel14_HouseholdIncomeFn.m

@@ -1,4 +1,4 @@
-function income=OLGModel14_HouseholdIncomeFn(h,sprime,s,z,e,agej,Jr,pension,r,w,P0,D,kappa_j,AccidentBeq,tau_l,tau_d,tau_cg,Lhscale)
+function income=OLGModel14_HouseholdIncomeFn(h,sprime,s,z,e,agej,Jr,pension,w,P0,D,kappa_j,AccidentBeq,r,tau_l,tau_d,tau_cg)
 % Replace assets with 'share holdings'
 % Get rid of progressive taxes
 % Add Lhnormalize
@@ -15,7 +15,7 @@
 if agej<Jr % If working age
     %consumption = labor income + accidental bequest + share holdings (including dividend) - capital gains tax - next period share holdings
     % income just is consumption but without subtracting the term for next period share holdings (-P*sprime)
-    income=(1-tau_l)*w*kappa_j*exp(z+e)*Lhscale*h+((1-tau_d)*D+P0)*(s+AccidentBeq) -tau_cg*(P0-Plag)*(s+AccidentBeq); 
+    income=(1-tau_l)*w*kappa_j*exp(z+e)*h+((1-tau_d)*D+P0)*(s+AccidentBeq) -tau_cg*(P0-Plag)*(s+AccidentBeq); 
 else % Retirement
     income=pension+((1-tau_d)*D+P0)*(s+AccidentBeq) -tau_cg*(P0-Plag)*(s+AccidentBeq);
 end

OLGModel14_HouseholdReturnFn.m

@@ -1,4 +1,4 @@
-function F=OLGModel14_HouseholdReturnFn(h,sprime,s,z,e,sigma,psi,eta,agej,Jr,J,pension,r,w,P0,D,kappa_j,warmglow1,warmglow2,AccidentBeq,tau_l,tau_d,tau_cg,Lhscale)
+function F=OLGModel14_HouseholdReturnFn(h,sprime,s,z,e,sigma,psi,eta,agej,Jr,J,pension,w,P0,D,kappa_j,warmglow1,warmglow2,AccidentBeq,r,tau_l,tau_d,tau_cg)
 % Replace assets with 'share holdings'
 % Get rid of progressive taxes
 % Add Lhnormalize
@@ -15,7 +15,7 @@
 F=-Inf;
 if agej<Jr % If working age
     %consumption = labor income + accidental bequest + share holdings (including dividend) - capital gains tax - next period share holdings
-    c=(1-tau_l)*w*kappa_j*exp(z+e)*Lhscale*h+((1-tau_d)*D+P0)*(s+AccidentBeq) -tau_cg*(P0-Plag)*(s+AccidentBeq)-P*sprime; 
+    c=(1-tau_l)*w*kappa_j*exp(z+e)*h+((1-tau_d)*D+P0)*(s+AccidentBeq) -tau_cg*(P0-Plag)*(s+AccidentBeq)-P*sprime; 
 else % Retirement
     c=pension+((1-tau_d)*D+P0)*(s+AccidentBeq) -tau_cg*(P0-Plag)*(s+AccidentBeq) - P*sprime;
 end