compute starting wealth as a function of CLI parameters and use it. f…

…ixes #275

sharkfin/utilities.py

@@ -1,10 +1,14 @@
 from HARK.core import distribute_params
 from HARK.distribution import Uniform
+from HARK.ConsumptionSaving.ConsPortfolioModel import SequentialPortfolioConsumerType
 import math
 import numpy as np
 import random
 from itertools import chain
 
+from scipy.interpolate import CubicSpline
+from scipy.optimize import fsolve
+
 
 # Distribution Utilities
 def update_return(dict1, dict2):
@@ -119,6 +123,12 @@ def combine_lognormal_rates(ror1, std1, ror2, std2):
     return ror3, sigma3
 
 
+def interp_func(x,y):
+    def func(z):
+        return np.interp(z, x, y)
+
+    return func
+
 ##### Lucas Pricing Equations
 
 import math
@@ -177,17 +187,108 @@ def lucas_expected_rate_of_return(pdr, dgr, dsd):
 
     return daily_ror, daily_std
 
-def expected_quarterly_returns(dgr, dst):
+def expected_quarterly_returns(DiscFac, CRRA, dgr, dst, days_per_quarter):
+    # DiscFac - quarterly discount factor
+    # CRRA - CRRA
     # dgr - daily dividend growth rate
     # dst - daily dividend standard deviation
 
     pdr = price_dividend_ratio_random_walk(
-        quarterly_params["DiscFac"],
-        annual_params["CRRA"],
+        DiscFac,
+        CRRA,
         dgr,
         dst,
-        60
+        days_per_quarter
     )
 
     (ror, sig) = lucas_expected_rate_of_return(pdr, dgr, dst)
-    return ror_quarterly(ror, 60), sig_quarterly(sig, 60)
+    return ror_quarterly(ror, days_per_quarter), sig_quarterly(sig, days_per_quarter)
+
+
+## Computing the target wealth
+
+def compute_target_wealth(
+    CRRA=6.0,
+    DiscFac=0.9,
+    RiskyAvg=1.08,
+    RiskyStd=0.20,
+    PermShkStd=[0.0],
+    PermGroFac=[1.0001],
+    UnempPrb=0.00
+):
+    agent_parameters = {}
+
+    agent_parameters["CRRA"] = CRRA
+    agent_parameters["DiscFac"] = DiscFac
+    agent_parameters["RiskyAvg"] = RiskyAvg
+    agent_parameters["RiskyStd"] = RiskyStd
+    agent_parameters["PermShkStd"] = PermShkStd
+    agent_parameters["PermGroFac"] = PermGroFac
+    agent_parameters["UnempPrb"] = UnempPrb
+    agent_parameters["LivPrb"] = [1.0]
+
+    agent = SequentialPortfolioConsumerType(**agent_parameters)
+
+    linear_roots, log_linear_roots, cubic_spline_roots = [], [], []
+
+    try:
+        agent.solve()
+        solved = True
+    except Exception as e:
+        solved = False
+
+        return solved, linear_roots, log_linear_roots, cubic_spline_roots
+
+    cFunc = agent.solution[0].cFuncAdj
+    ShareFunc = agent.solution[0].ShareFuncAdj
+
+    def expected_increase(ShareFunc, cFunc, mNrm):
+        share = ShareFunc(mNrm)
+        aNrm = mNrm - cFunc(mNrm)
+
+        mNrm_next = (
+            aNrm
+            * (
+                share * agent.parameters["RiskyAvg"]
+                + (1 - share) * agent.parameters["Rfree"]
+            )
+            + 1
+        )
+
+        gain = mNrm_next - aNrm
+        return gain
+
+    def expected_m_next(mNrm):
+        share = ShareFunc(mNrm)
+        aNrm = mNrm - cFunc(mNrm)
+        mNrm_next = (
+            aNrm
+            * (
+                share * agent.parameters["RiskyAvg"]
+                + (1 - share) * agent.parameters["Rfree"]
+            )
+            + 1
+        )
+
+        return mNrm_next
+
+    mNrm = np.linspace(0, 5, 1000)
+
+    # plt.plot(mNrm, cFunc(mNrm), label="c")
+
+    #plt.plot(mNrm, mNrm - expected_m_next(mNrm), label="m - E[m']")
+
+    linear_roots = fsolve(interp_func(mNrm, mNrm - expected_m_next(mNrm)), [mNrm[0]])
+    log_linear_roots = np.log(fsolve(interp_func(mNrm, mNrm - expected_m_next(mNrm)), [mNrm[0]]))
+    cubic_spline_roots = CubicSpline(mNrm,  mNrm - expected_m_next(mNrm)).roots()
+    print(f"m - E[m] linear interp roots: {linear_roots}")
+    print(f"m - E[m] log roots: {log_linear_roots}")
+    print(f"m - E[m] CubicSpine roots: {cubic_spline_roots}")
+
+    #plt.plot(mNrm, np.zeros_like(mNrm), label="0")
+
+    #plt.plot(mNrm, (mNrm - cFunc(mNrm)) * ShareFunc(mNrm), label ="wealth-into-market" )
+
+    #plt.legend()
+
+    return solved, linear_roots, log_linear_roots, cubic_spline_roots

macro/macro_parameters.py → simulate/macro_parameters.py

@@ -43,68 +43,3 @@
 ""
 quarterly_params["RiskyStd"]
 
-""
-quarterly_params
-
-###############################################################################
-# --- Computing the risky expectations from the dividend statistics
-
-from sharkfin.utilities import expected_quarterly_returns
-
-
-""
-
-
-
-""
-expected_quarterly_returns(1.0002, 0.011)
-
-""
-1.013848  ** 4
-
-""
-n = 100
-
-dgr_grid = np.linspace(1.000, 1.0005, n)
-dst_grid = np.linspace(0.007, 0.015, n)
-
-rors = np.zeros((n,n))
-sigs = np.zeros((n,n))
-
-for i in range(n):
-    for j in range(n):
-        try:
-            ror, sig = expected_quarterly_returns(dgr_grid[i], dst_grid[j])
-            rors[i,j] = ror - 0.0123
-            sigs[i,j] = sig - 0.1
-        except:
-            rors[i,j] = np.nan
-            sigs[i,j] = np.nan
-
-
-
-""
-rors[:20, :]
-
-""
-import matplotlib.pyplot as plt
-
-""
-error = np.absolute(ror) + np.absolute(sigs)
-
-""
-plt.imshow(error)
-plt.colorbar()
-
-###############################################################################
-# Looking for '0' on both of these plots below. Notice how close that puts us to the boundary.
-
-plt.imshow(rors)
-plt.colorbar()
-
-""
-plt.imshow(sigs)
-plt.colorbar()
-
-""
-

simulate/run_any_simulation.py

@@ -26,8 +26,14 @@
 from sharkfin.markets.ammps import ClientRPCMarket
 from sharkfin.population import SharkPopulation
 from sharkfin.simulation import AttentionSimulation, CalibrationSimulation
-from sharkfin.utilities import price_dividend_ratio_random_walk
+from sharkfin.utilities import (
+    price_dividend_ratio_random_walk,
+    lucas_expected_rate_of_return,
+    expected_quarterly_returns,
+    compute_target_wealth
+)
 
+from macro_parameters import quarterly_params
 
 class NpEncoder(json.JSONEncoder):
     def default(self, obj):
@@ -103,7 +109,7 @@ def default(self, obj):
 parser.add_argument(
     "--pop_aNrmInitMean",
     help="Log of initial mean asset levels for LUCAS0 population.",
-    default="0.5",
+    default=None,
 )
 
 parser.add_argument(
@@ -242,6 +248,34 @@ def run_chum_simulation(
 def env_param(name, default):
     return os.environ[name] if name in os.environ else default
 
+def target_log_wealth(
+        pop_CRRA,
+        pop_DiscFac,
+        dividend_growth_rate,
+        dividend_std,
+        days_per_quarter
+):
+
+    ror, sig = expected_quarterly_returns(
+        pop_DiscFac,
+        pop_CRRA,
+        dividend_growth_rate,
+        dividend_std,
+        days_per_quarter
+        )
+
+    solved, linear_roots, log_linear_roots, cubic_spline_roots = compute_target_wealth(
+        CRRA=pop_CRRA,
+        DiscFac=pop_DiscFac,
+        RiskyAvg=ror,
+        RiskyStd=sig,
+        PermShkStd=quarterly_params["PermShkStd"],
+        PermGroFac=quarterly_params["PermGroFac"],
+        UnempPrb=quarterly_params["UnempPrb"]
+    )
+
+    return np.log(linear_roots)[0]
+
 
 if __name__ == "__main__":
     args = parser.parse_args()
@@ -259,7 +293,7 @@ def env_param(name, default):
     population_name = str(args.population)
     pop_CRRA = float(args.pop_CRRA)
     pop_DiscFac = float(args.pop_DiscFac)
-    pop_aNrmInitMean = float(args.pop_aNrmInitMean)
+    pop_aNrmInitMean = float(args.pop_aNrmInitMean) if args.pop_aNrmInitMean is not None else None
 
     expectations_class_name = str(args.expectations)
     dividend_growth_rate = float(args.dividend_growth_rate)
@@ -269,6 +303,10 @@ def env_param(name, default):
     attention = float(args.attention)
     dphm = int(args.dphm)
 
+    if pop_aNrmInitMean is None:
+        pop_aNrmInitMean = target_log_wealth(pop_CRRA, pop_DiscFac, dividend_growth_rate, dividend_std, days_per_quarter)
+        print(f"Computed target wealth: {pop_aNrmInitMean}")
+
     # Memory-based FinanceModel arguments
     p1 = float(args.p1)
     p2 = float(args.p2)
@@ -330,13 +368,15 @@ def env_param(name, default):
     # random number generator with seed
     rng = np.random.default_rng(seed)
 
+    pdr = price_dividend_ratio_random_walk(
+            pop_DiscFac, pop_CRRA, dividend_growth_rate, dividend_std, days_per_quarter
+        )
+
     market_args = {
         "dividend_growth_rate": dividend_growth_rate,
         "dividend_std": dividend_std,
         "rng": rng,
-        "price_to_dividend_ratio": price_dividend_ratio_random_walk(
-            pop_DiscFac, pop_CRRA, dividend_growth_rate, dividend_std, days_per_quarter
-        ),
+        "price_to_dividend_ratio": pdr,
     }
 
     market_class = None
@@ -437,6 +477,9 @@ def env_param(name, default):
 
     with open(f"{filename}_sim_stats.txt", "w+") as f:
         sim_stats["filename"] = filename
+        sim_stats["dividend_std"] = dividend_std
+        sim_stats["pop_aNrmInitMean"] = pop_aNrmInitMean
+        sim_stats["price_dividend_ratio"] = pdr
         f.write(json.dumps(sim_stats, cls=NpEncoder))
 
     try: