{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# 税率一定の資産のランダム分配の帰結\n", "\n", "* Author: 黒木 玄\n", "* Date: 2019-03-20\n", "\n", "統計力学の教科書を見ると, 人口が多いとき, 資産を定額徴収してランダムに分配し続けると, 資産分布は指数分布に収束するという事実が解説されている.\n", "\n", "それでは, 資産を**定額**徴収(**税額**一定)するのではなく, **定率**徴収(**税率**一定)してランダム分配し続けるとどうなるだろうか?\n", "\n", "下の方のシミュレーションの結果を見れば, 税率一定のランダム分配を続けると, 資産分布はガンマ分布に収束し, 税率が**低い**ほど, 収束先の資産分布は**平等**に近付くことがわかる. この計算結果を初めて知ると直観に反する感じがして驚く." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" }, "toc": true }, "source": [ "

目次

\n", "
" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## シミュレーションとプロットのコード" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "dispimg (generic function with 1 method)" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Random: seed!\n", "\n", "using Base64\n", "dispimg(mime, file) = open(file) do f\n", " base64 = base64encode(f)\n", " display(\"text/html\", \"\"\"\"\"\")\n", "end" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### シミュレーションのコード" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "using Distributions\n", "\n", "function asset_update!(asset, taxfunc, niters)\n", " n = lastindex(asset)\n", " ur = UnitRange{Int32}(1:n)\n", " for iter in 1:niters\n", " k = rand(ur)\n", " l = rand(ur)\n", " d = taxfunc(asset[k])\n", " asset[k] -= d\n", " asset[l] += d\n", " end\n", "end\n", "\n", "function sim!(asset, taxfunc; niters=5*10^5, nframes=100, thin=niters÷nframes)\n", " n = lastindex(asset)\n", " T = eltype(asset)\n", " niters = thin * nframes\n", " Asset = zeros(T, n, nframes+1)\n", " Asset[:,1] .= asset\n", " for t in 1:nframes\n", " asset_update!(asset, taxfunc, thin)\n", " Asset[:,t+1] .= asset\n", " end\n", " Asset\n", "end\n", "\n", "function Elog(Asset)\n", " n, L = size(Asset)\n", " ElogA = Array{eltype(Asset),1}(undef, L)\n", " for l in 1:L\n", " ElogA[l] = mean(log(Asset[k,l]) for k in 1:n)\n", " end\n", " ElogA\n", "end\n", "\n", "struct ConstantRateTax{T<:AbstractFloat}; rate::T; end\n", "ConstantRateTax(; rate = 0.2) = ConstantRateTax(rate)\n", "(f::ConstantRateTax)(x) = f.rate * x" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### プロット用のコード" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "makegif (generic function with 1 method)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Plots\n", "gr()\n", "ENV[\"PLOTS_TEST\"] = \"true\"\n", "\n", "Plots.reset_defaults()\n", "gr(\n", " legend=false,\n", " titlefont=font(\"sans-serif\", 12),\n", " legendfont=font(\"sans-serif\", 8),\n", " guidefont=font(\"sans-serif\", 10),\n", " tickfont=font(\"sans-serif\", 8),\n", " markerstrokewidth=0,\n", ")\n", "\n", "function plot_Asset(t, Asset, ElogA; xmax=4.0, ymax=2.0, plottitle=\"\")\n", " n, L = size(Asset)\n", " fitdist = fit(Gamma, @view(Asset[:,t+1]))\n", " x = range(0, xmax, length=201)\n", "\n", " p1 = histogram(@view(Asset[:,t+1]), normed=true, alpha=0.5, linealpha=0, label=\"asset dist\")\n", " plot!(p1, xlim=(0,xmax), ylim=(0,ymax))\n", " plot!(p1, x, pdf.(fitdist, x), label=\"Gamma dist\", legend=true)\n", " if plottitle == \"\"\n", " title!(p1, \"t = $t\")\n", " else\n", " title!(p1, \"$plottitle (t = $t)\")\n", " end\n", "\n", " ElogAmin, ElogAmax = extrema(ElogA)\n", " ElogAmargin = 0.05*(ElogAmax - ElogAmin)\n", " ElogAmin -= ElogAmargin\n", " ElogAmax += ElogAmargin\n", " p2 = plot(0:t, ElogA[1:t+1])\n", " plot!(p2, [L], [0], xlim=(0,L-1), ylim=(ElogAmin, ElogAmax))\n", " title!(p2, \"mean of log(assets)\")\n", "\n", " plot(p1, p2, size=(700, 250))\n", "end\n", "\n", "function makegif(giffile, Asset; xmax=4.0, ymax=2.0, plottitle=\"\")\n", " N = size(Asset)[2]-1\n", " ElogA = Elog(Asset)\n", " \n", " # t = 0～L について1フレーム分の画像を作成\n", " update(t) = plot_Asset(t, Asset, ElogA; xmax=xmax, ymax=ymax, plottitle=plottitle)\n", " \n", " # t の動かす範囲の指定\n", " frames = [0;0;0;0;0;0;0;0;0:N;N;N;N;N;N;N;N;N;N:-1:0]\n", "\n", " # アニメーションオブジェクトの作成\n", " myanim = Animation()\n", "\n", " # アニメーションのすべてのフレームを作成\n", " for t in frames\n", " update(t)\n", " frame(myanim)\n", " end\n", " \n", " # GIFファイルの作成\n", " gif(myanim, giffile, fps=10)\n", "end" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## 税率一定の場合のシミュレーションの実行" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "### 資産税率90%" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0.387606 seconds (707.62 k allocations: 55.665 MiB, 8.53% gc time)\n" ] }, { "data": { "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ " 32.548658 seconds (96.35 M allocations: 3.813 GiB, 6.57% gc time)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "┌ Info: Saved animation to \n", "│ fn = C:\\Users\\genkuroki\\OneDrive\\msfd28\\images\\asset_dist_rate0.9.gif\n", "└ @ Plots C:\\Users\\genkuroki\\.julia\\packages\\Plots\\QYETN\\src\\animation.jl:90\n" ] } ], "source": [ "seed!(2018)\n", "n = 2^16\n", "asset = rand(Gamma(100, 0.01), n)\n", "rate = 0.9\n", "taxfunc = ConstantRateTax(rate=rate)\n", "@time Asset = sim!(asset, taxfunc, niters=10^6, nframes=40)\n", "\n", "giffile = \"images/asset_dist_rate$rate.gif\"\n", "@time makegif(giffile, Asset; ymax=10.0, plottitle=\"rate=$rate\")\n", "dispimg(\"image/gif\", giffile)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### 資産税率50%" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0.072987 seconds (57 allocations: 20.502 MiB, 23.77% gc time)\n" ] }, { "data": { "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ " 6.145228 seconds (40.12 M allocations: 1.197 GiB, 4.92% gc time)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "┌ Info: Saved animation to \n", "│ fn = C:\\Users\\genkuroki\\OneDrive\\msfd28\\images\\asset_dist_rate0.5.gif\n", "└ @ Plots C:\\Users\\genkuroki\\.julia\\packages\\Plots\\QYETN\\src\\animation.jl:90\n" ] } ], "source": [ "seed!(2018)\n", "n = 2^16\n", "asset = rand(Gamma(100, 0.01), n)\n", "rate = 0.5\n", "taxfunc = ConstantRateTax(rate=rate)\n", "@time Asset = sim!(asset, taxfunc, niters=10^6, nframes=40)\n", "\n", "giffile = \"images/asset_dist_rate$rate.gif\"\n", "@time makegif(giffile, Asset; ymax=4.0, plottitle=\"rate=$rate\")\n", "dispimg(\"image/gif\", giffile)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### 資産税率20%" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0.071026 seconds (1.98 k allocations: 20.615 MiB, 3.37% gc time)\n" ] }, { "data": { "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ " 6.539030 seconds (43.79 M allocations: 1.280 GiB, 3.02% gc time)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "┌ Info: Saved animation to \n", "│ fn = C:\\Users\\genkuroki\\OneDrive\\msfd28\\images\\asset_dist_rate0.2.gif\n", "└ @ Plots C:\\Users\\genkuroki\\.julia\\packages\\Plots\\QYETN\\src\\animation.jl:90\n" ] } ], "source": [ "seed!(2018)\n", "n = 2^16\n", "asset = rand(Exponential(1.0), n)\n", "rate = 0.2\n", "taxfunc = ConstantRateTax(rate=rate)\n", "@time Asset = sim!(asset, taxfunc, niters=15*10^5, nframes=40)\n", "\n", "giffile = \"images/asset_dist_rate$rate.gif\"\n", "@time makegif(giffile, Asset; ymax=1.0, plottitle=\"rate=$rate\")\n", "dispimg(\"image/gif\", giffile)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### 資産税率5%" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0.128894 seconds (58 allocations: 20.502 MiB, 2.74% gc time)\n" ] }, { "data": { "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ " 6.204482 seconds (42.38 M allocations: 1.249 GiB, 4.68% gc time)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "┌ Info: Saved animation to \n", "│ fn = C:\\Users\\genkuroki\\OneDrive\\msfd28\\images\\asset_dist_rate0.05.gif\n", "└ @ Plots C:\\Users\\genkuroki\\.julia\\packages\\Plots\\QYETN\\src\\animation.jl:90\n" ] } ], "source": [ "seed!(2018)\n", "n = 2^16\n", "asset = rand(Exponential(1.0), n)\n", "rate = 0.05\n", "taxfunc = ConstantRateTax(rate=rate)\n", "@time Asset = sim!(asset, taxfunc, niters=4*10^6, nframes=40)\n", "\n", "giffile = \"images/asset_dist_rate$rate.gif\"\n", "@time makegif(giffile, Asset; ymax=2.0, plottitle=\"rate=$rate\")\n", "dispimg(\"image/gif\", giffile)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### 資産税率1%" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0.479512 seconds (1.98 k allocations: 20.615 MiB, 1.00% gc time)\n" ] }, { "data": { "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ " 6.454085 seconds (40.10 M allocations: 1.196 GiB, 3.01% gc time)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "┌ Info: Saved animation to \n", "│ fn = C:\\Users\\genkuroki\\OneDrive\\msfd28\\images\\asset_dist_rate0.01.gif\n", "└ @ Plots C:\\Users\\genkuroki\\.julia\\packages\\Plots\\QYETN\\src\\animation.jl:90\n" ] } ], "source": [ "seed!(2018)\n", "n = 2^16\n", "asset = rand(Exponential(1.0), n)\n", "rate = 0.01\n", "taxfunc = ConstantRateTax(rate=rate)\n", "@time Asset = sim!(asset, taxfunc, niters=2*10^7, nframes=40)\n", "\n", "giffile = \"images/asset_dist_rate$rate.gif\"\n", "@time makegif(giffile, Asset; ymax=4.0, plottitle=\"rate=$rate\")\n", "dispimg(\"image/gif\", giffile)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 税額一定の場合のシミュレーションの実行" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "税額が一定だと分布の収束先は指数分布になる.\n", "\n", "こちらについては統計力学の教科書などによく書かれている." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "struct ConstantTax{T<:AbstractFloat}; tax::T; end\n", "ConstantTax(; tax = 0.05) = ConstantTax(tax)\n", "(f::ConstantTax)(x) = ifelse(f.tax ≥ x, zero(x), f.tax)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 2.682780 seconds (254.28 k allocations: 32.727 MiB, 0.20% gc time)\n" ] }, { "data": { "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ " 6.135301 seconds (39.71 M allocations: 1.182 GiB, 4.57% gc time)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "┌ Info: Saved animation to \n", "│ fn = C:\\Users\\genkuroki\\OneDrive\\msfd28\\images\\asset_dist_tax0.05.gif\n", "└ @ Plots C:\\Users\\genkuroki\\.julia\\packages\\Plots\\QYETN\\src\\animation.jl:90\n" ] } ], "source": [ "seed!(2018)\n", "n = 2^16\n", "asset = rand(Gamma(100, 0.01), n)\n", "tax = 0.05\n", "taxfunc = ConstantTax(tax=tax)\n", "@time Asset = sim!(asset, taxfunc, niters=10^8, nframes=40)\n", "\n", "giffile = \"images/asset_dist_tax$tax.gif\"\n", "@time makegif(giffile, Asset; ymax=4.0, plottitle=\"tax=$tax\")\n", "dispimg(\"image/gif\", giffile)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "_draft": { "nbviewer_url": "https://gist.github.com/82a163e0204beb41bb2bd9ae91b0483e" }, "gist": { "data": { "description": "msfd28/税率一定の資産のランダム分配の帰結.ipynb", "public": true }, "id": "82a163e0204beb41bb2bd9ae91b0483e" }, "kernelspec": { "display_name": "Julia 1.1.0", "language": "julia", "name": "julia-1.1" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.1.0" }, "toc": { "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": true, "title_cell": "目次", "title_sidebar": "目次", "toc_cell": true, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }