{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Symbolic computation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Pkg.add(\"SymPy\")\n",
    "using SymPy, LinearAlgebra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(a, b, c)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@vars a b c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(x, y, z)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#we can optionally specify domain with e.g.\n",
    "x,y,z = symbols(\"x, y, z\", real=true) # also integer, positive, negative..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Algebra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}x + 2 y\\end{equation*}"
      ],
      "text/plain": [
       "x + 2⋅y"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr = x+2y # note that we do not need to specify 2*y.. this is Julia specific"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}x \\left(x + 2 y\\right)\\end{equation*}"
      ],
      "text/plain": [
       "x⋅(x + 2⋅y)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr2 = x*expr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}x^{2} + 2 x y\\end{equation*}"
      ],
      "text/plain": [
       " 2        \n",
       "x  + 2⋅x⋅y"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expanded_expr = SymPy.expand(x*expr)  # Attenction that without the star it gives wrong results!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}x \\left(x + 2 y\\right)\\end{equation*}"
      ],
      "text/plain": [
       "x⋅(x + 2⋅y)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "factor(expanded_expr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}\\left(x + 1\\right)^{2}\\end{equation*}"
      ],
      "text/plain": [
       "       2\n",
       "(x + 1) "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = (x+1)^2 # power is ^ and not ** .. again, Julia specific"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}x^{2} - 2 x + 1\\end{equation*}"
      ],
      "text/plain": [
       " 2          \n",
       "x  - 2⋅x + 1"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = x^2-2x+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}- x^{2} + 2 x + \\left(x + 1\\right)^{2} - 1\\end{equation*}"
      ],
      "text/plain": [
       "   2                2    \n",
       "- x  + 2⋅x + (x + 1)  - 1"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f-g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}4 x\\end{equation*}"
      ],
      "text/plain": [
       "4⋅x"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simplify(f-g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}2 z^{3} + 1\\end{equation*}"
      ],
      "text/plain": [
       "   3    \n",
       "2⋅z  + 1"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr3 = subs(expr, x=>1, y=>z^3) # substitution can be either numeric or symbolic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[ \\left[ \\begin{array}{r}- \\frac{x}{2}\\end{array} \\right] \\]"
      ],
      "text/plain": [
       "1-element Array{Sym,1}:\n",
       " -x/2"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve(expr2,y) # by default solve the equation f(.) = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[ \\left[ \\begin{array}{r}0\\\\- 2 y\\end{array} \\right] \\]"
      ],
      "text/plain": [
       "2-element Array{Sym,1}:\n",
       "    0\n",
       " -2*y"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve(expr2,x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}x \\left(x + 2 y\\right) = 2\\end{equation*}"
      ],
      "text/plain": [
       "x⋅(x + 2⋅y) = 2"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Eq(expr2,2) # create an equation object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[ \\left[ \\begin{array}{r}- \\frac{x}{2} + \\frac{1}{x}\\end{array} \\right] \\]"
      ],
      "text/plain": [
       "1-element Array{Sym,1}:\n",
       " -x/2 + 1/x"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve(Eq(expr2,2),y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1-element Array{Tuple{Sym,Sym},1}:\n",
       " (1, -1/2)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve([expr,g],[x,y]) # solve systems of equations in multiple variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix algebra\n",
    " - matrices are just Julian matrices with symbolic entries\n",
    " - constructing matrices then follows Julia's conventions\n",
    " - compared with \"plain\" Julia\" here we can do symbolic operations, not just numeric ones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\left[ \\begin{array}{rr}1&4\\\\2&5\\\\x&x\\end{array}\\right]\\]"
      ],
      "text/plain": [
       "3×2 Array{Sym,2}:\n",
       " 1  4\n",
       " 2  5\n",
       " x  x"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = [[1,2,x] [4,5,x]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\left[ \\begin{array}{rrr}x&3&5\\\\2&4&6\\end{array}\\right]\\]"
      ],
      "text/plain": [
       "2×3 Array{Sym,2}:\n",
       " x  3  5\n",
       " 2  4  6"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = [[x,2] [3,4] [5,6]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\left[ \\begin{array}{rrr}x + 8&19&29\\\\2 x + 10&26&40\\\\x^{2} + 2 x&7 x&11 x\\end{array}\\right]\\]"
      ],
      "text/plain": [
       "3×3 Array{Sym,2}:\n",
       "     x + 8   19    29\n",
       "  2*x + 10   26    40\n",
       " x^2 + 2*x  7*x  11*x"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M*N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}6 x^{2} + 11 x \\left(26 x + 208\\right) - 11 x \\left(38 x + 190\\right) - 7 x \\left(40 x + 320\\right) + 7 x \\left(58 x + 290\\right) + 12 x\\end{equation*}"
      ],
      "text/plain": [
       "   2                                                                          \n",
       "6⋅x  + 11⋅x⋅(26⋅x + 208) - 11⋅x⋅(38⋅x + 190) - 7⋅x⋅(40⋅x + 320) + 7⋅x⋅(58⋅x + \n",
       "\n",
       "           \n",
       "290) + 12⋅x"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "det(M*N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[ \\left[ \\begin{array}{r}x\\\\y\\\\2\\end{array} \\right] \\]"
      ],
      "text/plain": [
       "3-element Array{Sym,1}:\n",
       " x\n",
       " y\n",
       " 2"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = [a,b,3]\n",
    "e = [x,y,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}x \\overline{a} + y \\overline{b} + 6\\end{equation*}"
      ],
      "text/plain": [
       "  _     _    \n",
       "x⋅a + y⋅b + 6"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k = d' * e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[\\left[ \\begin{array}{rrr}a x&a y&2 a\\\\b x&b y&2 b\\\\3 x&3 y&6\\end{array}\\right]\\]"
      ],
      "text/plain": [
       "3×3 Array{Sym,2}:\n",
       " a*x  a*y  2*a\n",
       " b*x  b*y  2*b\n",
       " 3*x  3*y    6"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "j = d * e'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}a x^{2} + 3 c - 4 x\\end{equation*}"
      ],
      "text/plain": [
       "   2            \n",
       "a⋅x  + 3⋅c - 4⋅x"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = a*x^2-4x+3c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}2 a x - 4\\end{equation*}"
      ],
      "text/plain": [
       "2⋅a⋅x - 4"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(f,x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[ \\left[ \\begin{array}{r}\\frac{2}{a}\\end{array} \\right] \\]"
      ],
      "text/plain": [
       "1-element Array{Sym,1}:\n",
       " 2/a"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solve(diff(f,x),x) # find FOC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}\\frac{a x^{3}}{3} + 3 c x - 2 x^{2}\\end{equation*}"
      ],
      "text/plain": [
       "   3               \n",
       "a⋅x               2\n",
       "──── + 3⋅c⋅x - 2⋅x \n",
       " 3                 "
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate(f,x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}\\frac{56 a}{3} + 6 c - 24\\end{equation*}"
      ],
      "text/plain": [
       "56⋅a           \n",
       "──── + 6⋅c - 24\n",
       " 3             "
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "integrate(f,(x,2,4)) # definite integral between x=2 and x=4 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## An application: derive Faustmann formula"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(t, k, r, p, NPV)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@vars t k r p NPV       # time, planting costs, interest rate, timber price, Net Present Value (as symbol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(S,)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@symfuns S              # generic function for stock"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Net Present Value \n",
    "- $PV ~=~ \\frac{b}{(1+r)^t-1}$\n",
    "- $b ~=~  {p S\\{t\\} - k (1+r)^t}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}\\frac{- k + p S{\\left (t \\right )} e^{- r t}}{1 - e^{- r t}}\\end{equation*}"
      ],
      "text/plain": [
       "             -r⋅t\n",
       "-k + p⋅S(t)⋅ℯ    \n",
       "─────────────────\n",
       "         -r⋅t    \n",
       "    1 - ℯ        "
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "NPVf = ( p * S(t) * exp(-r*t)-k)/ (1-exp(-r*t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}- \\frac{r \\left(- k + p S{\\left (t \\right )} e^{- r t}\\right) e^{- r t}}{\\left(1 - e^{- r t}\\right)^{2}} + \\frac{- p r S{\\left (t \\right )} e^{- r t} + p e^{- r t} \\frac{d}{d t} S{\\left (t \\right )}}{1 - e^{- r t}}\\end{equation*}"
      ],
      "text/plain": [
       "                                            -r⋅t      -r⋅t d       \n",
       "    ⎛             -r⋅t⎞  -r⋅t   - p⋅r⋅S(t)⋅ℯ     + p⋅ℯ    ⋅──(S(t))\n",
       "  r⋅⎝-k + p⋅S(t)⋅ℯ    ⎠⋅ℯ                                  dt      \n",
       "- ─────────────────────────── + ───────────────────────────────────\n",
       "                     2                            -r⋅t             \n",
       "          ⎛     -r⋅t⎞                        1 - ℯ                 \n",
       "          ⎝1 - ℯ    ⎠                                              "
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Derivative of NPV with respect to time\n",
    "NPVt = diff(NPVf,t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}- \\frac{r \\left(- k + p S{\\left (t \\right )} e^{- r t}\\right) e^{- r t}}{\\left(1 - e^{- r t}\\right)^{2}} + \\frac{- p r S{\\left (t \\right )} e^{- r t} + p e^{- r t} \\frac{d}{d t} S{\\left (t \\right )}}{1 - e^{- r t}} = 0\\end{equation*}"
      ],
      "text/plain": [
       "                                            -r⋅t      -r⋅t d           \n",
       "    ⎛             -r⋅t⎞  -r⋅t   - p⋅r⋅S(t)⋅ℯ     + p⋅ℯ    ⋅──(S(t))    \n",
       "  r⋅⎝-k + p⋅S(t)⋅ℯ    ⎠⋅ℯ                                  dt          \n",
       "- ─────────────────────────── + ─────────────────────────────────── = 0\n",
       "                     2                            -r⋅t                 \n",
       "          ⎛     -r⋅t⎞                        1 - ℯ                     \n",
       "          ⎝1 - ℯ    ⎠                                                  "
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Equation of derivative of NPV equal to zero\n",
    "NPVteq = Eq(NPVt,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\[ \\left[ \\begin{array}{r}- \\frac{r \\left(k - p S{\\left (t \\right )}\\right) e^{r t}}{e^{r t} - 1}\\end{array} \\right] \\]"
      ],
      "text/plain": [
       "1-element Array{Sym,1}:\n",
       " -r*(k - p*S(t))*exp(r*t)/(exp(r*t) - 1)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Rewriting the equation with the value derivative as dependent variable  (V' = pS')\n",
    "Vt = solve(NPVteq,p * sympy.Derivative(S(t), t)) # V' = f(pS,r,k,t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}\\frac{- k + p S{\\left (t \\right )} e^{- r t}}{1 - e^{- r t}} = NPV\\end{equation*}"
      ],
      "text/plain": [
       "             -r⋅t      \n",
       "-k + p⋅S(t)⋅ℯ          \n",
       "───────────────── = NPV\n",
       "         -r⋅t          \n",
       "    1 - ℯ              "
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Setting the \"symbol\" NPV equal to the NPV function\n",
    "# (because we want to keep in the final output the result in term of NPV)\n",
    "NPVeq = Eq(NPVf,NPV)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}\\left(- NPV e^{r t} + NPV + p S{\\left (t \\right )}\\right) e^{- r t}\\end{equation*}"
      ],
      "text/plain": [
       "⎛       r⋅t               ⎞  -r⋅t\n",
       "⎝- NPV⋅ℯ    + NPV + p⋅S(t)⎠⋅ℯ    "
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solve the equation NPV==NPVb with respect to k, which has a unique solution as a function of NPV\n",
    "k_NPV = solve(NPVeq, k)[1] # k = f(NPV)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}- \\frac{r \\left(- p S{\\left (t \\right )} + \\left(- NPV e^{r t} + NPV + p S{\\left (t \\right )}\\right) e^{- r t}\\right) e^{r t}}{e^{r t} - 1}\\end{equation*}"
      ],
      "text/plain": [
       "   ⎛          ⎛       r⋅t               ⎞  -r⋅t⎞  r⋅t \n",
       "-r⋅⎝-p⋅S(t) + ⎝- NPV⋅ℯ    + NPV + p⋅S(t)⎠⋅ℯ    ⎠⋅ℯ    \n",
       "──────────────────────────────────────────────────────\n",
       "                        r⋅t                           \n",
       "                       ℯ    - 1                       "
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Now substitute k\n",
    "VtNPV = subs(Vt[1], k=>k_NPV) # V' = f(NPV)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}r \\left(NPV + p S{\\left (t \\right )}\\right)\\end{equation*}"
      ],
      "text/plain": [
       "r⋅(NPV + p⋅S(t))"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Now simplify\n",
    "VtNPV2 = simplify(VtNPV) # pS' = f(pS, NPV, r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{equation*}NPV r + p r S{\\left (t \\right )}\\end{equation*}"
      ],
      "text/plain": [
       "NPV⋅r + p⋅r⋅S(t)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Now expand\n",
    "SymPy.expand(VtNPV2) # pS' = f(pS, NPV, r), the Faustmann formula"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.2.0",
   "language": "julia",
   "name": "julia-1.2"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}