{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "closed-glenn", "metadata": { "tags": [] }, "outputs": [], "source": [ "import torch\n", "from torch.autograd.functional import jacobian\n", "import itertools" ] }, { "cell_type": "markdown", "id": "naval-ivory", "metadata": {}, "source": [ "# Setup Functions\n", "## General CompositionFunctions" ] }, { "cell_type": "code", "execution_count": 2, "id": "italian-enforcement", "metadata": {}, "outputs": [], "source": [ "### Set up functions to compose functions \n", "# These functions will \n", "# - compose two functions together\n", "# - compose a function to itself n times.\n", "\n", "def compose(f,g):\n", " return lambda *args: f(g(*args))\n", "\n", "def compose_recursive_functions(fn,n):\n", " #takes a function fn and composes it with itself n times\n", " #Returns each of the iterations of the functions.\n", " \n", " #Set base conditions\n", " out_func = None\n", " out_func_list =[]\n", "\n", " #build the compositions of functions\n", " for f in itertools.repeat(fn, n):\n", " #if first iteration\n", " if out_func == None:\n", " out_func = f\n", " else:\n", " out_func = compose(f,out_func)\n", "\n", " #append the ou\n", " out_func_list.append(out_func)\n", " \n", " return out_func_list" ] }, { "cell_type": "markdown", "id": "fancy-tucson", "metadata": {}, "source": [ "## Setup functions related to the problem" ] }, { "cell_type": "code", "execution_count": 15, "id": "outside-arrangement", "metadata": {}, "outputs": [], "source": [ "### Background functions\n", "\n", "\n", "def survival(stock, debris):\n", " #Gompertz distribution for simplicity\n", " #commonly used with saturation\n", " #TODO: ACTUALLY DERIVE A SURVIVAL FUNCTION. THIS IS JUST A PLACEHOLDER. PROBABLY SHOULD BE AN EXPONENTIAL DISTRIBUTION\n", " \n", " #eta = 1.0/(SCALING@stock)\n", " #b = 1/debris\n", " #return 1 - ( b*eta*torch.exp(eta+b*stock-eta*torch.exp(b*stock)))\n", " return 1 - torch.exp(-SCALING * stock-debris)\n", "\n", "def test_launch(stock, debris):\n", " return torch.ones(5, requires_grad=True)\n", "\n", "def laws_of_motion(stock, debris, launches):\n", " \n", " new_stock = stock*survival(stock,debris) + launches#(stock,debris) #TODO: Launches will become a function (neural network)\n", " \n", " #TODO: Currently Ignoring autocatalysis\n", " new_debris = (1-DELTA)*debris + LAUNCH_DEBRIS_RATE * launches.sum() + COLLISION_DEBRIS_RATE*(1-survival(stock,debris)) @ stock\n", " \n", " return (new_stock, new_debris)\n", "\n", "#This is not a good specification of the profit function, but it will work for now.\n", "def profit(x):\n", " return UTIL_WEIGHTS @ x" ] }, { "cell_type": "code", "execution_count": 16, "id": "romance-generation", "metadata": {}, "outputs": [], "source": [ "def single_transition(item_to_iterate, laws_motion, profit,launch, stocks, debris ):\n", " #TODO: change launch from a direct tensor, to a function.\n", " \n", " #Calculate the inverse\n", " bA = BETA * jacobian(laws_motion, (stocks,debris,launch))[0][0]\n", " #TODO: figure out some diagnostics for this section\n", " \n", " \n", " return bA.inverse() @ (item_to_iterate - jacobian(profit,stocks))\n", "\n", "# This function wraps the single transition and handles updating dates etc.\n", "def transition_wrapper(data_in):\n", " #unpack states and functions\n", " stocks, debris,profit, launch, laws_motion, item_to_transition = data_in\n", " \n", " #Calculate new states\n", " new_stocks, new_debris = laws_motion(stocks,debris,launch)\n", " \n", " #WARNING: RECURSION: You may break your head...\n", " #This gets the transition of the value function derivatives over time.\n", " transitioned = single_transition(\n", " item_to_transition, #item to iterate, i.e. the derivatives of the value function\n", " #functions\n", " laws_motion, \n", " profit, \n", " launch, #TODO: reimplement with launch as a function\n", " stocks, debris #states\n", " )\n", " \n", " #collects the data back together for return, including the updated state variables\n", " data_out = new_stocks, new_debris, profit, launch, laws_motion, transitioned\n", " \n", " return data_out" ] }, { "cell_type": "markdown", "id": "fluid-parks", "metadata": {}, "source": [ "# Actual calculations" ] }, { "cell_type": "code", "execution_count": 17, "id": "changing-january", "metadata": {}, "outputs": [], "source": [ "#set states\n", "stocks = torch.ones(5)\n", "#Last one is different\n", "stocks[-1] = 0.5\n", "#now add the tracking requirement in place\n", "stocks.requires_grad_()\n", "\n", "#Setup Debris\n", "debris = torch.tensor([2.2],requires_grad=True)\n", "\n", "#CHANGE LATER: Launch is currently a value, should be a function (i.e. neural network)\n", "launch = torch.ones(5, requires_grad=True)\n", "\n", "#compose the functions together.\n", "base_data = (stocks,debris, profit, launch, laws_of_motion, torch.ones(5, requires_grad=True))\n", "\n", "#Parameters\n", "SCALING = torch.ones(5)\n", "DELTA = 0.9 \n", "LAUNCH_DEBRIS_RATE = 0.005\n", "COLLISION_DEBRIS_RATE = 0.0007\n", "UTIL_WEIGHTS = torch.tensor([1,-0.2,0,0,0])\n", "BETA = 0.95" ] }, { "cell_type": "code", "execution_count": 19, "id": "dominant-boost", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "tensor([0.0000, 1.2632, 1.0526, 1.0526, 1.0892], grad_fn=)\n", "tensor([-0.9519, 1.3928, 1.0020, 1.0020, 1.0575], grad_fn=)\n", "tensor([-1.8565, 1.5150, 0.9530, 0.9530, 0.9882], grad_fn=)\n", "tensor([-2.8103, 1.6872, 0.9376, 0.9376, 0.9474], grad_fn=)\n", "tensor([-3.8626, 1.9131, 0.9505, 0.9505, 0.9408], grad_fn=)\n", "tensor([-5.0235, 2.1830, 0.9819, 0.9819, 0.9598], grad_fn=)\n", "tensor([-6.2860, 2.4869, 1.0247, 1.0247, 0.9951], grad_fn=)\n", "tensor([-7.6403, 2.8175, 1.0746, 1.0746, 1.0403], grad_fn=)\n", "tensor([-9.0802, 3.1712, 1.1293, 1.1293, 1.0918], grad_fn=)\n", "tensor([-10.6035, 3.5462, 1.1879, 1.1879, 1.1477],\n", " grad_fn=)\n", "tensor([-12.2108, 3.9422, 1.2501, 1.2501, 1.2075],\n", " grad_fn=)\n", "tensor([-13.9045, 4.3597, 1.3157, 1.3157, 1.2708],\n", " grad_fn=)\n", "tensor([-15.6882, 4.7995, 1.3849, 1.3849, 1.3375],\n", " grad_fn=)\n", "tensor([-17.5662, 5.2625, 1.4577, 1.4577, 1.4079],\n", " grad_fn=)\n", "tensor([-19.5432, 5.7500, 1.5345, 1.5345, 1.4820],\n", " grad_fn=)\n", "tensor([-21.6244, 6.2631, 1.6152, 1.6152, 1.5600],\n", " grad_fn=)\n", "tensor([-23.8151, 6.8033, 1.7002, 1.7002, 1.6421],\n", " grad_fn=)\n", "tensor([-26.1212, 7.3719, 1.7897, 1.7897, 1.7285],\n", " grad_fn=)\n", "tensor([-28.5486, 7.9704, 1.8839, 1.8839, 1.8195],\n", " grad_fn=)\n", "tensor([-31.1038, 8.6004, 1.9831, 1.9831, 1.9152],\n", " grad_fn=)\n", "tensor([-33.7934, 9.2636, 2.0874, 2.0874, 2.0160],\n", " grad_fn=)\n", "tensor([-36.6247, 9.9617, 2.1973, 2.1973, 2.1221],\n", " grad_fn=)\n", "tensor([-39.6049, 10.6965, 2.3129, 2.3129, 2.2338],\n", " grad_fn=)\n", "tensor([-42.7420, 11.4700, 2.4347, 2.4347, 2.3514],\n", " grad_fn=)\n", "tensor([-46.0442, 12.2842, 2.5628, 2.5628, 2.4751],\n", " grad_fn=)\n", "tensor([-49.5202, 13.1413, 2.6977, 2.6977, 2.6054],\n", " grad_fn=)\n", "tensor([-53.1792, 14.0435, 2.8397, 2.8397, 2.7425],\n", " grad_fn=)\n", "tensor([-57.0307, 14.9931, 2.9891, 2.9891, 2.8869],\n", " grad_fn=)\n", "tensor([-61.0850, 15.9928, 3.1465, 3.1465, 3.0388],\n", " grad_fn=)\n", "tensor([-65.3526, 17.0450, 3.3121, 3.3121, 3.1988],\n", " grad_fn=)\n", "tensor([-69.8449, 18.1526, 3.4864, 3.4864, 3.3671],\n", " grad_fn=)\n", "tensor([-74.5736, 19.3186, 3.6699, 3.6699, 3.5443],\n", " grad_fn=)\n", "tensor([-79.5511, 20.5459, 3.8630, 3.8630, 3.7309],\n", " grad_fn=)\n", "tensor([-84.7907, 21.8378, 4.0664, 4.0664, 3.9272],\n", " grad_fn=)\n", "tensor([-90.3060, 23.1976, 4.2804, 4.2804, 4.1339],\n", " grad_fn=)\n", "tensor([-96.1115, 24.6291, 4.5057, 4.5057, 4.3515],\n", " grad_fn=)\n", "tensor([-102.2227, 26.1359, 4.7428, 4.7428, 4.5805],\n", " grad_fn=)\n", "tensor([-108.6555, 27.7220, 4.9924, 4.9924, 4.8216],\n", " grad_fn=)\n", "tensor([-115.4268, 29.3916, 5.2552, 5.2552, 5.0754],\n", " grad_fn=)\n", "tensor([-122.5545, 31.1490, 5.5318, 5.5318, 5.3425],\n", " grad_fn=)\n", "tensor([-130.0574, 32.9990, 5.8229, 5.8229, 5.6237],\n", " grad_fn=)\n", "tensor([-137.9552, 34.9463, 6.1294, 6.1294, 5.9197],\n", " grad_fn=)\n", "tensor([-146.2686, 36.9961, 6.4520, 6.4520, 6.2313],\n", " grad_fn=)\n", "tensor([-155.0196, 39.1538, 6.7916, 6.7916, 6.5592],\n", " grad_fn=)\n", "tensor([-164.2312, 41.4251, 7.1490, 7.1490, 6.9044],\n", " grad_fn=)\n", "tensor([-173.9275, 43.8158, 7.5253, 7.5253, 7.2678],\n", " grad_fn=)\n", "tensor([-184.1343, 46.3325, 7.9213, 7.9213, 7.6503],\n", " grad_fn=)\n", "tensor([-194.8782, 48.9815, 8.3382, 8.3382, 8.0530],\n", " grad_fn=)\n", "tensor([-206.1876, 51.7701, 8.7771, 8.7771, 8.4768],\n", " grad_fn=)\n", "tensor([-218.0922, 54.7053, 9.2391, 9.2391, 8.9230],\n", " grad_fn=)\n" ] } ], "source": [ "#calculate results for first 5 iterations\n", "for f in compose_recursive_functions(transition_wrapper,50):\n", " result = f(base_data)\n", " print(result[5])" ] }, { "cell_type": "markdown", "id": "unnecessary-architect", "metadata": {}, "source": [ "Note how this fails on the last few iterations.\n", "I need to get better model functions (profit, laws_of_motion, etc) together to test this out.\n", "Alternatively, I can check for areas where the determinant of $A$ is zero, possibly by doing some sort of grid search?\n", "\n", "Maybe with a standard RBC model?\n", "\n", "Also, maybe I can create a `Model` class that upon construction will capture the necesary constants, functions, etc." ] }, { "cell_type": "code", "execution_count": null, "id": "varying-organization", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.8" } }, "nbformat": 4, "nbformat_minor": 5 }