Plotly javascript loaded.
\n", "To load again call
init_notebook(true)\n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "Plots.PlotlyJSBackend()" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plotlyjs()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define a custom pretty printing function:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "prettyprint (generic function with 1 method)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "prettyprint(arr) = Base.showarray(STDOUT, arr, false, header=false)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 1: Reducible chain " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider the Markov chain given by the following stochastic matrix, taken from Exercise 3 in Jarvis and Shier (where the actual values of non-zero probabilities are not important):" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6×6 Array{Float64,2}:\n", " 1.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 1.0 0.0\n", " 0.0 0.0 0.333333 0.333333 0.333333 0.0\n", " 0.5 0.0 0.0 0.0 0.0 0.5\n", " 0.0 0.5 0.0 0.0 0.5 0.0\n", " 0.5 0.0 0.0 0.5 0.0 0.0" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "P = zeros(6, 6)\n", "P[1, 1] = 1\n", "P[2, 5] = 1\n", "P[3, 3], P[3, 4], P[3, 5] = 1/3, 1/3, 1/3\n", "P[4, 1], P[4, 6] = 1/2, 1/2\n", "P[5, 2], P[5, 5] = 1/2, 1/2\n", "P[6, 1], P[6, 4] = 1/2, 1/2\n", "P" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create a MarkovChain instance:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Discrete Markov Chain\n", "stochastic matrix of type Array{Float64,2}:\n", "[1.0 0.0 … 0.0 0.0; 0.0 0.0 … 1.0 0.0; … ; 0.0 0.5 … 0.5 0.0; 0.5 0.0 … 0.0 0.0]" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mc1 = MarkovChain(P)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Classification of states" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This Markov chain is reducible:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "false" ], "text/plain": [ "false" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "is_irreducible(mc1)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "4" ], "text/plain": [ "4" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "length(communication_classes(mc1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Determine the communication classes:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4-element Array{Array{Int64,1},1}:\n", " [1] \n", " [2, 5]\n", " [4, 6]\n", " [3] " ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "communication_classes(mc1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Classify the states of this Markov chain:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Array{Array{Int64,1},1}:\n", " [1] \n", " [2, 5]" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "recurrent_classes(mc1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Obtain a list of the recurrent states:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Int64,1}:\n", " 1\n", " 2\n", " 5" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "recurrent_states = vcat(recurrent_classes(mc1)...)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Obtain a list of the transient states:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Int64,1}:\n", " 3\n", " 4\n", " 6" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "transient_states = setdiff(collect(1:n_states(mc1)), recurrent_states)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A Markov chain is reducible (i.e., its directed graph is not strongly connected)\n", "if and only if by symmetric permulations of rows and columns,\n", "its transition probability matrix is written in the form (\"canonical form\")\n", "$$\n", "\\begin{pmatrix}\n", "U & 0 \\\\\n", "W & V\n", "\\end{pmatrix},\n", "$$\n", "where $U$ and $W$ are square matrices.\n", "\n", "Such a form for `mc1` is obtained by the following:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6×6 Array{Float64,2}:\n", " 1.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 1.0 0.0 0.0 0.0\n", " 0.0 0.5 0.5 0.0 0.0 0.0\n", " 0.0 0.0 0.333333 0.333333 0.333333 0.0\n", " 0.5 0.0 0.0 0.0 0.0 0.5\n", " 0.5 0.0 0.0 0.0 0.5 0.0" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "permutation = vcat(recurrent_states, transient_states)\n", "mc1.p[permutation, permutation]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This Markov chain is aperiodic\n", "(i.e., the least common multiple of the periods of the recurrent sub-chains is one):" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "true" ], "text/plain": [ "true" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "is_aperiodic(mc1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Indeed, each of the sub-chains corresponding to the recurrent classes has period $1$,\n", "i.e., every recurrent state is aperiodic:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Period of the sub-chain\n", " 1.0\n", " = 1\n", "Period of the sub-chain\n", " 0.0 1.0\n", " 0.5 0.5\n", " = 1\n" ] } ], "source": [ "for recurrent_class in recurrent_classes(mc1)\n", " sub_matrix = mc1.p[recurrent_class, recurrent_class]\n", " d = period(MarkovChain(sub_matrix))\n", " println(\"Period of the sub-chain\")\n", " prettyprint(sub_matrix)\n", " println(\"\\n = $d\")\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Stationary distributions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For each recurrent class $C$, there is a unique stationary distribution $\\psi^C$\n", "such that $\\psi^C_i > 0$ for all $i \\in C$ and $\\psi^C_i = 0$ otherwise.\n", "`MarkovChain.stationary_distributions` returns\n", "these unique stationary distributions for the recurrent classes.\n", "Any stationary distribution is written as a convex combination of these distributions." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Array{Array{Float64,1},1}:\n", " [1.0, 0.0, 0.0, 0.0, 0.0, 0.0] \n", " [0.0, 0.333333, 0.0, 0.0, 0.666667, 0.0]" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "stationary_distributions(mc1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These are indeed stationary distributions:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/plain": [ "1×6 RowVector{Float64,Array{Float64,1}}:\n", " 1.0 0.0 0.0 0.0 0.0 0.0" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "stationary_distributions(mc1)[1]'* mc1.p" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1×6 RowVector{Float64,Array{Float64,1}}:\n", " 0.0 0.333333 0.0 0.0 0.666667 0.0" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "stationary_distributions(mc1)[2]'* mc1.p" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot these distributions." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "draw_histogram" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", "Plot the given distribution.\n", "\"\"\"\n", "function draw_histogram(distribution;\n", " title=\"\", xlabel=\"\", ylabel=\"\", ylim=(0, 1),\n", " show_legend=false, show_grid=false)\n", " \n", " n = length(distribution)\n", " p = bar(collect(1:n),\n", " distribution,\n", " xlim=(0.5, n+0.5),\n", " ylim=ylim,\n", " title=title,\n", " xlabel=xlabel,\n", " ylabel=ylabel,\n", " xticks=1:n,\n", " legend=show_legend,\n", " grid=show_grid)\n", " \n", " return p\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Stationary distribution for the recurrent class:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "titles = [\"$recurrent_class\"\n", " for recurrent_class in recurrent_classes(mc1)]\n", "\n", "ps = []\n", "for (title, dist) in zip(titles, stationary_distributions(mc1))\n", " push!(ps, draw_histogram(dist, title=title, xlabel=\"States\"))\n", "end\n", "\n", "plot(ps..., layout=(1, 2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simulation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us simulate our Markov chain `mc1`.\n", "The `simualte` method generates a sample path\n", "of length given by the first argument, `ts_length`,\n", "with an initial state as specified by an optional argument `init`;\n", "if not specified, the initial state is randomly drawn. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A sample path from state `1`:\n", "\n", "(Note: Transposing the output array here is just for visualization.)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/plain": [ "1×50 RowVector{Int64,Array{Int64,1}}:\n", " 1 1 1 1 1 1 1 1 1 1 1 1 1 … 1 1 1 1 1 1 1 1 1 1 1 1" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "simulate(mc1, 50, init=1)'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As is clear from the transition matrix `P`,\n", "if it starts at state `1`, the chain stays there forever,\n", "i.e., `1` is an absorbing state, a state that constitutes a singleton recurrent class." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Start with state `2`:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1×50 RowVector{Int64,Array{Int64,1}}:\n", " 2 5 2 5 5 5 2 5 2 5 5 2 5 … 5 2 5 2 5 5 2 5 2 5 2 5" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "simulate(mc1, 50, init=2)'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can observe that the chain stays in the recurrent class $\\{2,5\\}$\n", "and visits states `2` and `5` with certain frequencies." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If `init` is not specified, the initial state is randomly chosen:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1×50 RowVector{Int64,Array{Int64,1}}:\n", " 2 5 2 5 2 5 5 2 5 5 5 5 2 … 5 2 5 5 2 5 2 5 5 5 5 5" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "simulate(mc1, 50)'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Time series averages" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let us compute the frequency distribution along a sample path, given by\n", "$$\n", "\\frac{1}{t} \\sum_{\\tau=0}^{t-1} \\mathbf{1}\\{X_{\\tau} = s\\}\n", "\\quad (s \\in S).\n", "$$" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "time_series_dist (generic function with 2 methods)" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", "Return the distribution of visits by a sample path of length t\n", "of mc with an initial state init.\n", "\"\"\"\n", "function time_series_dist(mc, ts::Array{Int}; init=rand(1:n_states(mc)))\n", " t_max = maximum(ts)\n", " ts_size = length(ts)\n", " \n", " X = simulate(mc, t_max, init=init)\n", " dists = Array{Float64}(n_states(mc), ts_size)\n", " bins = 1:n_states(mc)+1\n", " for (i, t) in enumerate(ts)\n", " h = fit(Histogram, X[1:t], bins, closed=:left)\n", " dists[:, i] = h.weights / t\n", " end\n", " \n", " return dists\n", "end\n", "\n", "function time_series_dist(mc, t::Int; init=rand(1:n_states(mc)))\n", " return time_series_dist(mc, [t], init=init)[:, 1]\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is a frequency distribution along a sample path, of length 100,\n", "from initial state `2`, which is a recurrent state:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6-element Array{Float64,1}:\n", " 0.0 \n", " 0.33\n", " 0.0 \n", " 0.0 \n", " 0.67\n", " 0.0 " ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "time_series_dist(mc1, 100, init=2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Length 10,000:" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6-element Array{Float64,1}:\n", " 0.0 \n", " 0.3311\n", " 0.0 \n", " 0.0 \n", " 0.6689\n", " 0.0 " ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "time_series_dist(mc1, 10^4, init=2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The distribution becomes close to the stationary distribution `(0, 1/3, 0, 0, 2/3, 0)`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the frequency distributions for a couple of different time lengths:" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "plot_time_series_dists (generic function with 1 method)" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function plot_time_series_dists(mc, ts; init=rand(1:n_states(mc)), layout=(1,length(ts)))\n", " dists = time_series_dist(mc, ts, init=init)\n", " titles = [\"t=$t\" for t in ts]\n", "\n", " ps = []\n", " for (i, title) in enumerate(titles)\n", " p = draw_histogram(dists[:, i], title=title, xlabel=\"States\")\n", " push!(ps, p)\n", " end\n", " plot(ps..., layout=layout)\n", "end" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "init = 2\n", "ts = [5, 10, 50, 100]\n", "plot_time_series_dists(mc1, ts, init=init)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Start with state `3`,\n", "which is a transient state:" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "init = 3\n", "ts = [5, 10, 50, 100]\n", "plot_time_series_dists(mc1, ts, init=init)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run the above cell several times;\n", "you will observe that the limit distribution differs across sample paths.\n", "Sometimes the state is absorbed into the recurrent class $\\{1\\}$,\n", "while other times it is absorbed into the recurrent class $\\{2,5\\}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some sample path with init=3" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "init = 3\n", "ts = [5, 10, 50, 100]\n", "\n", "plot_time_series_dists(mc1, ts, init=init)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another sample path with init=3" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot_time_series_dists(mc1, ts, init=init)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In fact,\n", "for almost every sample path of a finite Markov chain $\\{X_t\\}$,\n", "for some recurrent class $C$ we have\n", "$$\n", "\\frac{1}{t} \\sum_{\\tau=0}^{t-1} \\mathbf{1}\\{X_{\\tau} = s\\} \\to \\psi^C[s]\n", "\\quad \\text{as $t \\to \\infty$}\n", "$$\n", "for all states $s$,\n", "where $\\psi^C$ is the stationary distribution associated with the recurrent class $C$.\n", "\n", "If the initial state $s_0$ is a recurrent state,\n", "then the recurrent class $C$ above is the one that contains $s_0$,\n", "while if it is a transient state,\n", "then the recurrent class to which the convergence occurs depends on the sample path." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us simulate with the remaining states, `4`, `5`, and `6`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Time series distributions for t=100" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "inits = [4, 5, 6]\n", "t = 100\n", "\n", "ps = []\n", "for init in inits\n", " p = draw_histogram(\n", " time_series_dist(mc1, t, init=init),\n", " title=\"Initial state = $init\",\n", " xlabel=\"States\")\n", " push!(ps, p)\n", "end\n", "\n", "plot(ps..., layout=(1, 3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Cross sectional averages" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, let us repeat the simulation for many times (say, 10,000 times)\n", "and obtain the distribution of visits to each state at a given time period `t`.\n", "That is, we want to simulate the marginal distribution at time `t`." ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": true }, "outputs": [], "source": [ "function QuantEcon.simulate(mc, ts_length, num_reps::Int=10^4; init=1)\n", " X = Array{Int}(ts_length, num_reps)\n", " for i in 1:num_reps\n", " X[:, i] = simulate(mc, ts_length, init=init)\n", " end\n", " return X\n", "end" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "cross_sectional_dist (generic function with 4 methods)" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", "Return the distribution of visits at time T by num_reps times of simulation\n", "of mc with an initial state init.\n", "\"\"\"\n", "function cross_sectional_dist(mc, ts::Array{Int}, num_reps=10^4; init=rand(1:n_states(mc)))\n", " t_max = maximum(ts)\n", " ts_size = length(ts)\n", " dists = Array{Float64}(n_states(mc), ts_size)\n", " X = simulate(mc, t_max+1, num_reps, init=init)\n", " bins = 1:n_states(mc)+1\n", " for (i, t) in enumerate(ts)\n", " h = fit(Histogram, X[t, :], bins, closed=:left)\n", " dists[:, i] = h.weights / num_reps\n", " end\n", " return dists\n", "end\n", "\n", "function cross_sectional_dist(mc, t::Int, num_reps=10^4; init=rand(1:n_states(mc)))\n", " return cross_sectional_dist(mc, [t], num_reps, init=init)[:, 1]\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Start with state `2`:" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6-element Array{Float64,1}:\n", " 0.0 \n", " 0.3396\n", " 0.0 \n", " 0.0 \n", " 0.6604\n", " 0.0 " ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "init = 2\n", "t = 10\n", "cross_sectional_dist(mc1, t, init=init)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6-element Array{Float64,1}:\n", " 0.0 \n", " 0.3398\n", " 0.0 \n", " 0.0 \n", " 0.6602\n", " 0.0 " ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "t = 100\n", "cross_sectional_dist(mc1, t, init=init)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The distribution is close to the stationary distribution (0, 1/3, 0, 0, 2/3, 0)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the simulated marginal distribution at `t` for some values of `t`." ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "plot_cross_sectional_dists (generic function with 2 methods)" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function plot_cross_sectional_dists(mc, ts, num_reps=10^4; init=1)\n", " dists = cross_sectional_dist(mc, ts, num_reps, init=init)\n", " titles = map(t -> \"t=$t\", ts)\n", "\n", " ps = []\n", " for (i, title) in enumerate(titles)\n", " p = draw_histogram(dists[:, i], title=title, xlabel=\"States\")\n", " push!(ps, p)\n", " end\n", "\n", " plot(ps..., layout=(1, length(ts)))\n", "end" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "init = 2\n", "ts = [2, 3, 5, 10]\n", "plot_cross_sectional_dists(mc1, ts, init=init)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Starting with a transient state `3`:" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6-element Array{Float64,1}:\n", " 0.4965\n", " 0.1658\n", " 0.0 \n", " 0.0027\n", " 0.333 \n", " 0.002 " ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "init = 3\n", "t = 10\n", "cross_sectional_dist(mc1, t, init=init)" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6-element Array{Float64,1}:\n", " 0.5044\n", " 0.1636\n", " 0.0 \n", " 0.0 \n", " 0.332 \n", " 0.0 " ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "t = 100\n", "dist = cross_sectional_dist(mc1, t, init=init)" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "draw_histogram(dist,\n", " title=\"Cross sectional distribution at t=$t with init=$init\",\n", " xlabel=\"States\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Observe that the distribution is close to a convex combination of\n", "the stationary distributions `(1, 0, 0, 0, 0, 0)` and `(0, 1/3, 0, 0, 2/3, 0)`,\n", "which is a stationary distribution itself." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How the simulated marginal distribution evolves:" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "init = 3\n", "ts = [2, 3, 5, 10]\n", "plot_cross_sectional_dists(mc1, ts, init=init)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since our Markov chain is aperiodic (i.e., every recurrent class is aperiodic),\n", "the marginal disribution at time $t$ converges as $t \\to \\infty$\n", "to some stationary distribution,\n", "and the limit distribution depends on the initial state,\n", "according to the probabilities that the state is absorbed into the recurrent classes." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For initial states `4`, `5`, and `6`:" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "inits = [4, 5, 6]\n", "t = 10\n", "\n", "ps = []\n", "for init in inits\n", " p = draw_histogram(cross_sectional_dist(mc1, t, init=init),\n", " title=\"Initial state = $init\",\n", " xlabel=\"States\")\n", " push!(ps, p)\n", "end\n", "\n", "plot(ps..., layout=(1, length(inits)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Powers of $P$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The marginal distributions at time $t$ are obtained by $P^t$." ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "P^10 =\n", " 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.333984 0.0 0.0 0.666016 0.0 \n", " 0.498072 0.166792 1.69351e-5 0.000767702 0.3332 0.00115155 \n", " 0.999023 0.0 0.0 0.000976563 0.0 0.0 \n", " 0.0 0.333008 0.0 0.0 0.666992 0.0 \n", " 0.999023 0.0 0.0 0.0 0.0 0.000976563\n", "P^20 =\n", " 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.333334 0.0 0.0 0.666666 0.0 \n", " 0.499998 0.166667 2.86797e-10 7.6271e-7 0.333333 1.14407e-6\n", " 0.999999 0.0 0.0 9.53674e-7 0.0 0.0 \n", " 0.0 0.333333 0.0 0.0 0.666667 0.0 \n", " 0.999999 0.0 0.0 0.0 0.0 9.53674e-7\n", "P^30 =\n", " 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.333333 0.0 0.0 0.666667 0.0 \n", " 0.5 0.166667 4.85694e-15 7.45054e-10 0.333333 1.11758e-9 \n", " 1.0 0.0 0.0 9.31323e-10 0.0 0.0 \n", " 0.0 0.333333 0.0 0.0 0.666667 0.0 \n", " 1.0 0.0 0.0 0.0 0.0 9.31323e-10\n" ] } ], "source": [ "ts = [10, 20, 30]\n", "for t in ts\n", " P_t = mc1.p^t\n", " println(\"P^$t =\")\n", " prettyprint(P_t)\n", " println()\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the canonical form:" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Q =\n", " 1.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 1.0 0.0 0.0 0.0\n", " 0.0 0.5 0.5 0.0 0.0 0.0\n", " 0.0 0.0 0.333333 0.333333 0.333333 0.0\n", " 0.5 0.0 0.0 0.0 0.0 0.5\n", " 0.5 0.0 0.0 0.0 0.5 0.0\n", "Q^10 =\n", " 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.333984 0.666016 0.0 0.0 0.0 \n", " 0.0 0.333008 0.666992 0.0 0.0 0.0 \n", " 0.498072 0.166792 0.3332 1.69351e-5 0.000767702 0.00115155 \n", " 0.999023 0.0 0.0 0.0 0.000976563 0.0 \n", " 0.999023 0.0 0.0 0.0 0.0 0.000976563\n", "Q^20 =\n", " 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.333334 0.666666 0.0 0.0 0.0 \n", " 0.0 0.333333 0.666667 0.0 0.0 0.0 \n", " 0.499998 0.166667 0.333333 2.86797e-10 7.6271e-7 1.14407e-6\n", " 0.999999 0.0 0.0 0.0 9.53674e-7 0.0 \n", " 0.999999 0.0 0.0 0.0 0.0 9.53674e-7\n", "Q^30 =\n", " 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.333333 0.666667 0.0 0.0 0.0 \n", " 0.0 0.333333 0.666667 0.0 0.0 0.0 \n", " 0.5 0.166667 0.333333 4.85694e-15 7.45054e-10 1.11758e-9 \n", " 1.0 0.0 0.0 0.0 9.31323e-10 0.0 \n", " 1.0 0.0 0.0 0.0 0.0 9.31323e-10" ] } ], "source": [ "Q = mc1.p[permutation, permutation]\n", "println(\"Q =\")\n", "prettyprint(Q)\n", "for t in ts\n", " Q_t = Q^t\n", " println(\"\\nQ^$t =\")\n", " prettyprint(Q_t)\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Observe that the first three rows, which correspond to the recurrent states,\n", "are close to the stationary distributions associated with the corresponding recurrent classes." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 2: Periodic chain" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider the Markov chain given by the following stochastic matrix,\n", "taken from Exercise 9 (see also Exercise 11) in Jarvis and Shier\n", "(where the actual values of non-zero probabilities are not important):" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "10×10 Array{Float64,2}:\n", " 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.5 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0\n", " 0.0 0.333333 0.333333 0.0 0.0 0.0 0.0 0.333333 0.0 0.0\n", " 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.5 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0\n", " 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "P = zeros(10, 10)\n", "P[1, 4] = 1\n", "P[2, [1, 5]] = 1/2\n", "P[3, 7] = 1\n", "P[4, [2, 3, 8]] = 1/3\n", "P[5, 4] = 1\n", "P[6, 5] = 1\n", "P[7, 4] = 1\n", "P[8, [7, 9]] = 1/2\n", "P[9, 10] = 1\n", "P[10, 6] = 1\n", "P" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Discrete Markov Chain\n", "stochastic matrix of type Array{Float64,2}:\n", "[0.0 0.0 … 0.0 0.0; 0.5 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 1.0; 0.0 0.0 … 0.0 0.0]" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mc2 = MarkovChain(P)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This Markov chain is irreducible:" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "data": { "text/html": [ "true" ], "text/plain": [ "true" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "is_irreducible(mc2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This Markov chain is periodic:" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "data": { "text/html": [ "false" ], "text/plain": [ "false" ] }, "execution_count": 48, "metadata": {}, "output_type": "execute_result" } ], "source": [ "is_aperiodic(mc2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Its period, which we denote by $d$:" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [ { "data": { "text/html": [ "3" ], "text/plain": [ "3" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d = period(mc2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Cyclic classes are:\n", "\n", "Note: The function to identify the cyclic classes is only implemented in Python version and not available in Julia version.\n", "The result below is retrieved from the paper introduced above: [Graph-Theoretic Analysis of Finite Markov Chains]((http://www.ces.clemson.edu/~shierd/Shier/markov.pdf)" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Array{Int64,1},1}:\n", " [1, 5, 7, 9]\n", " [4, 10] \n", " [2, 3, 6, 8]" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cyclic_classes = [[1, 5, 7, 9], [4, 10], [2, 3, 6, 8]]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Cyclic normal form" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If a Markov chain is periodic with period $d \\geq 2$,\n", "then its transition probability matrix is written in the form (\"cyclic normal form\")\n", "$$\n", "\\begin{pmatrix}\n", "0 & P_1 & 0 & 0 & \\cdots & 0 \\\\\n", "0 & 0 & P_2 & 0 & \\cdots & 0 \\\\\n", "0 & 0 & 0 & P_3 & \\cdots & 0 \\\\\n", "\\vdots & \\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n", "0 & 0 & 0 & 0 & \\cdots & P_{d-1} \\\\\n", "P_d & 0 & 0 & 0 & \\cdots & 0\n", "\\end{pmatrix}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Represent our Markov chain in cyclic normal form:" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "10×10 Array{Float64,2}:\n", " 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.333333 0.333333 0.0 0.333333\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 \n", " 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 " ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "permutation = vcat(cyclic_classes...)\n", "Q = mc2.p[permutation, permutation]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Re-define the Markov chain with the above matrix `Q`:" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Discrete Markov Chain\n", "stochastic matrix of type Array{Float64,2}:\n", "[0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 1.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]" ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mc2 = MarkovChain(Q)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Obtain the block components $P_1, \\cdots, P_{d}$:" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Array{Int64,1},1}:\n", " [1, 2, 3, 4] \n", " [5, 6] \n", " [7, 8, 9, 10]" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cyclic_classes = [[1, 2, 3, 4], [5, 6], [7, 8, 9, 10]]" ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "P_1 =\n", " 1.0 0.0\n", " 1.0 0.0\n", " 1.0 0.0\n", " 0.0 1.0\n", "P_2 =\n", " 0.333333 0.333333 0.0 0.333333\n", " 0.0 0.0 1.0 0.0 \n", "P_3 =\n", " 0.5 0.5 0.0 0.0\n", " 0.0 0.0 1.0 0.0\n", " 0.0 1.0 0.0 0.0\n", " 0.0 0.0 0.5 0.5\n" ] } ], "source": [ "P_blocks = []\n", "\n", "for i in 1:d\n", " push!(P_blocks, mc2.p[cyclic_classes[(i-1)%d+1], :][:, cyclic_classes[i%d+1]])\n", " println(\"P_$i =\")\n", " prettyprint(P_blocks[i])\n", " println()\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$P^d$ is block diagonal:" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0.166667 0.166667 0.5 0.166667 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.166667 0.166667 0.5 0.166667 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.166667 0.166667 0.5 0.166667 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.833333 0.166667 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.333333 0.333333 0.0 0.333333\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.333333 0.333333 0.0 0.333333\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.333333 0.333333 0.0 0.333333\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.166667 0.166667 0.5 0.166667" ] } ], "source": [ "P_power_d = mc2.p^d\n", "prettyprint(P_power_d)" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1st diagonal block of P^d =\n", " 0.166667 0.166667 0.5 0.166667\n", " 0.166667 0.166667 0.5 0.166667\n", " 0.166667 0.166667 0.5 0.166667\n", " 0.0 1.0 0.0 0.0 \n", "2nd diagonal block of P^d =\n", " 0.833333 0.166667\n", " 1.0 0.0 \n", "3rd diagonal block of P^d =\n", " 0.333333 0.333333 0.0 0.333333\n", " 0.333333 0.333333 0.0 0.333333\n", " 0.333333 0.333333 0.0 0.333333\n", " 0.166667 0.166667 0.5 0.166667\n" ] } ], "source": [ "P_power_d_blocks = []\n", "ordinals = [\"1st\", \"2nd\", \"3rd\"]\n", "\n", "for (i, ordinal) in enumerate(ordinals)\n", " push!(P_power_d_blocks, P_power_d[cyclic_classes[i], :][:, cyclic_classes[i]])\n", " println(\"$ordinal diagonal block of P^d =\")\n", " prettyprint(P_power_d_blocks[i])\n", " println()\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The $i$th diagonal block of $P^d$ equals $P_i P_{i+1} \\cdots P_d P_1 \\cdots P_{i-1}$:" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "P_1 P_2 P_3 =\n", " 0.166667 0.166667 0.5 0.166667\n", " 0.166667 0.166667 0.5 0.166667\n", " 0.166667 0.166667 0.5 0.166667\n", " 0.0 1.0 0.0 0.0 \n", "P_2 P_3 P_1 =\n", " 0.833333 0.166667\n", " 1.0 0.0 \n", "P_3 P_1 P_2 =\n", " 0.333333 0.333333 0.0 0.333333\n", " 0.333333 0.333333 0.0 0.333333\n", " 0.333333 0.333333 0.0 0.333333\n", " 0.166667 0.166667 0.5 0.166667\n" ] } ], "source": [ "products = []\n", "\n", "for i in 1:d\n", " R = P_blocks[i]\n", " string = \"P_$i\"\n", " for j in 1:d-1\n", " R = R * P_blocks[(i+j-1)%d+1]\n", " string *= \" P_$((i+j-1)%d+1)\"\n", " end\n", " push!(products, R)\n", " println(string, \" =\")\n", " prettyprint(R)\n", " println()\n", "end" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "true\n", "true\n", "true\n" ] } ], "source": [ "for (matrix0, matrix1) in zip(P_power_d_blocks, products)\n", " println(matrix0 == matrix1)\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Stationary distributions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Markov chain `mc2` has a unique stationary distribution,\n", "which we denote by $\\psi$:" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [ { "data": { "text/html": [ "1" ], "text/plain": [ "1" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" } ], "source": [ "length(stationary_distributions(mc2))" ] }, { "cell_type": "code", "execution_count": 60, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "10-element Array{Float64,1}:\n", " 0.047619 \n", " 0.0952381\n", " 0.142857 \n", " 0.047619 \n", " 0.285714 \n", " 0.047619 \n", " 0.0952381\n", " 0.0952381\n", " 0.047619 \n", " 0.0952381" ] }, "execution_count": 60, "metadata": {}, "output_type": "execute_result" } ], "source": [ "psi = stationary_distributions(mc2)[1]" ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ "draw_histogram(psi,\n", " title=\"Stationary distribution\",\n", " xlabel=\"States\",\n", " ylim=(0, 0.35))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Obtain the stationary distributions $\\psi^1, \\ldots, \\psi^{d}$\n", "each associated with the diagonal blocks of $P^d$:" ] }, { "cell_type": "code", "execution_count": 62, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "psi^1 =\n", "[0.142857, 0.285714, 0.428571, 0.142857]\n", "psi^2 =\n", "[0.857143, 0.142857]\n", "psi^3 =\n", "[0.285714, 0.285714, 0.142857, 0.285714]\n" ] } ], "source": [ "psi_s = []\n", "\n", "for i in 1:d\n", " push!(psi_s, stationary_distributions(MarkovChain(P_power_d_blocks[i]))[1])\n", " println(\"psi^$i =\")\n", " println(psi_s[i])\n", "end" ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 63, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ps = []\n", "\n", "for i in 1:d\n", " psi_i_full_dim = zeros(n_states(mc2))\n", " psi_i_full_dim[cyclic_classes[i]] = psi_s[i]\n", " p = draw_histogram(psi_i_full_dim,\n", " title=\"ψ^$i\",\n", " xlabel=\"States\")\n", " push!(ps, p)\n", "end\n", "\n", "plot(ps..., layout=(1, 3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Verify that $\\psi^{i+1} = \\psi^i P_i$:" ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "psi^1 P_1 =\n", "[0.857143 0.142857]\n", "psi^2 P_2 =\n", "[0.285714 0.285714 0.142857 0.285714]\n", "psi^3 P_3 =\n", "[0.142857 0.285714 0.428571 0.142857]\n" ] } ], "source": [ "for i in 1:d\n", " println(\"psi^$i P_$i =\")\n", " println(psi_s[i]' * P_blocks[i])\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Verify that $\\psi = (\\psi^1 + \\cdots + \\psi^d)/d$:" ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "10-element Array{Float64,1}:\n", " 0.047619 \n", " 0.0952381\n", " 0.142857 \n", " 0.047619 \n", " 0.285714 \n", " 0.047619 \n", " 0.0952381\n", " 0.0952381\n", " 0.047619 \n", " 0.0952381" ] }, "execution_count": 65, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Right hand side of the above identity\n", "rhs = zeros(n_states(mc2))\n", "\n", "for i in 1:d\n", " rhs[cyclic_classes[i]] = psi_s[i]\n", "end\n", "rhs /= d\n", "rhs" ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [ { "data": { "text/html": [ "2.7755575615628914e-17" ], "text/plain": [ "2.7755575615628914e-17" ] }, "execution_count": 66, "metadata": {}, "output_type": "execute_result" } ], "source": [ "maximum(abs.(psi - rhs))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "meaning right hand side is very close to $\\psi$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Powers of $P$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since the Markov chain in consideration is periodic,\n", "the marginal distribution does not converge, but changes periodically." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us compute the powers of the transition probability matrix (in cyclic normal form):" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Print $P^1, P^2, \\ldots, P^d$:" ] }, { "cell_type": "code", "execution_count": 67, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "P^1 =\n", " 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.333333 0.333333 0.0 0.333333\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 \n", " 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 \n", "P^2 =\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.333333 0.333333 0.0 0.333333\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.333333 0.333333 0.0 0.333333\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.333333 0.333333 0.0 0.333333\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 \n", " 0.166667 0.166667 0.5 0.166667 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.5 0.5 0.0 0.0 0.0 0.0 \n", "P^3 =\n", " 0.166667 0.166667 0.5 0.166667 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.166667 0.166667 0.5 0.166667 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.166667 0.166667 0.5 0.166667 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.833333 0.166667 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.333333 0.333333 0.0 0.333333\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.333333 0.333333 0.0 0.333333\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.333333 0.333333 0.0 0.333333\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.166667 0.166667 0.5 0.166667\n", "P^4 =\n", " 0.0 0.0 0.0 0.0 0.833333 0.166667 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.833333 0.166667 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.833333 0.166667 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.277778 0.277778 0.166667 0.277778\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.333333 0.333333 0.0 0.333333\n", " 0.166667 0.166667 0.5 0.166667 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.166667 0.166667 0.5 0.166667 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.166667 0.166667 0.5 0.166667 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0833333 0.583333 0.25 0.0833333 0.0 0.0 0.0 0.0 0.0 0.0 \n" ] } ], "source": [ "for i in 1:d+1\n", " println(\"P^$i =\")\n", " prettyprint(mc2.p^i)\n", " println()\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Print $P^{2d}$, $P^{4d}$, and $P^{6d}$:" ] }, { "cell_type": "code", "execution_count": 68, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "P^6 =\n", " 0.138889 0.305556 0.416667 0.138889 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.138889 0.305556 0.416667 0.138889 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.138889 0.305556 0.416667 0.138889 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.166667 0.166667 0.5 0.166667 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.861111 0.138889 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.833333 0.166667 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.277778 0.277778 0.166667 0.277778\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.277778 0.277778 0.166667 0.277778\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.277778 0.277778 0.166667 0.277778\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.305556 0.305556 0.0833333 0.305556\n", "P^12 =\n", " 0.142747 0.286265 0.428241 0.142747 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.142747 0.286265 0.428241 0.142747 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.142747 0.286265 0.428241 0.142747 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.143519 0.282407 0.430556 0.143519 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.857253 0.142747 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.856481 0.143519 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.285494 0.285494 0.143519 0.285494\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.285494 0.285494 0.143519 0.285494\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.285494 0.285494 0.143519 0.285494\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.286265 0.286265 0.141204 0.286265\n", "P^18 =\n", " 0.142854 0.28573 0.428562 0.142854 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.142854 0.28573 0.428562 0.142854 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.142854 0.28573 0.428562 0.142854 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.142876 0.285622 0.428627 0.142876 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.857146 0.142854 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.857124 0.142876 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.285708 0.285708 0.142876 0.285708\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.285708 0.285708 0.142876 0.285708\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.285708 0.285708 0.142876 0.285708\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.28573 0.28573 0.142811 0.28573 \n" ] } ], "source": [ "for i in [k*d for k in [2, 4, 6]]\n", " println(\"P^$i =\")\n", " prettyprint(mc2.p^i)\n", " println()\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$P^{kd}$ converges as $k \\to \\infty$ to a matrix that contains $\\psi^1, \\ldots, \\psi^d$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Print $P^{kd+1}, \\ldots, P^{kd+d}$ with $k = 10$ for example:" ] }, { "cell_type": "code", "execution_count": 69, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "P^31 =\n", " 0.0 0.0 0.0 0.0 0.857143 0.142857 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.857143 0.142857 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.857143 0.142857 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.857143 0.142857 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.285714 0.285714 0.142857 0.285714\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.285714 0.285714 0.142857 0.285714\n", " 0.142857 0.285714 0.428571 0.142857 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.142857 0.285714 0.428571 0.142857 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.142857 0.285714 0.428571 0.142857 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.142857 0.285714 0.428571 0.142857 0.0 0.0 0.0 0.0 0.0 0.0 \n", "P^32 =\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.285714 0.285714 0.142857 0.285714\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.285714 0.285714 0.142857 0.285714\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.285714 0.285714 0.142857 0.285714\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.285714 0.285714 0.142857 0.285714\n", " 0.142857 0.285714 0.428571 0.142857 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.142857 0.285714 0.428571 0.142857 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.857143 0.142857 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.857143 0.142857 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.857143 0.142857 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.857143 0.142857 0.0 0.0 0.0 0.0 \n", "P^33 =\n", " 0.142857 0.285714 0.428571 0.142857 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.142857 0.285714 0.428571 0.142857 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.142857 0.285714 0.428571 0.142857 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.142857 0.285714 0.428571 0.142857 0.0 0.0 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.857143 0.142857 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.857143 0.142857 0.0 0.0 0.0 0.0 \n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.285714 0.285714 0.142857 0.285714\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.285714 0.285714 0.142857 0.285714\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.285714 0.285714 0.142857 0.285714\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.285714 0.285714 0.142857 0.285714\n" ] } ], "source": [ "for i in 10*d+1:11*d\n", " println(\"P^$i =\")\n", " prettyprint(mc2.p^i)\n", " println()\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But $P^i$ itself does not converge." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the frequency distribution of visits to the states\n", "along a sample path starting at state `1`:" ] }, { "cell_type": "code", "execution_count": 70, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "10-element Array{Float64,1}:\n", " 0.0449\n", " 0.1008\n", " 0.1372\n", " 0.0505\n", " 0.2828\n", " 0.0505\n", " 0.0951\n", " 0.0908\n", " 0.0505\n", " 0.0969" ] }, "execution_count": 70, "metadata": {}, "output_type": "execute_result" } ], "source": [ "init = 1\n", "dist = time_series_dist(mc2, 10^4, init=init)" ] }, { "cell_type": "code", "execution_count": 71, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 71, "metadata": {}, "output_type": "execute_result" } ], "source": [ "draw_histogram(dist,\n", " title=\"Time series distribution with init=$init\",\n", " xlabel=\"States\", ylim=(0, 0.35))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Observe that the distribution is close to the (unique) stationary distribution $\\psi$." ] }, { "cell_type": "code", "execution_count": 72, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "10-element Array{Float64,1}:\n", " 0.047619 \n", " 0.0952381\n", " 0.142857 \n", " 0.047619 \n", " 0.285714 \n", " 0.047619 \n", " 0.0952381\n", " 0.0952381\n", " 0.047619 \n", " 0.0952381" ] }, "execution_count": 72, "metadata": {}, "output_type": "execute_result" } ], "source": [ "psi" ] }, { "cell_type": "code", "execution_count": 73, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 73, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bar(psi, legend=false, xlabel=\"States\", title=\"ψ\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, plot the simulated marginal distributions\n", "at $T = 10d+1, \\ldots, 11d, 11d+1, \\ldots, 12d$ with initial state `1`:" ] }, { "cell_type": "code", "execution_count": 74, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ] }, "execution_count": 74, "metadata": {}, "output_type": "execute_result" } ], "source": [ "init = 1\n", "k = 10\n", "ts = [k*d + i + 1 for i in 1:2*d]\n", "num_reps = 10^2\n", "dists = cross_sectional_dist(mc2, ts, num_reps, init=init)\n", "\n", "ps = []\n", "for (i, t) in enumerate(ts)\n", " p = draw_histogram(dists[:, i],\n", " title=\"t = $t\")\n", " push!(ps, p)\n", "end\n", "plot(ps..., layout=(2, d))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare these with the rows of $P^{10d+1}, \\ldots, P^{10d+d}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 3: Nearly completely decomposable chain" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider the Markov chain given by the following stochastic matrix $P^{\\varepsilon}$,\n", "parameterized by $\\varepsilon$:" ] }, { "cell_type": "code", "execution_count": 75, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "P_epsilon (generic function with 2 methods)" ] }, "execution_count": 75, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function P_epsilon(eps, p=0.5)\n", " P = [1-(p+eps) p eps;\n", " p 1-(p+eps) eps;\n", " eps eps 1-2*eps]\n", " return P\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If $\\varepsilon = 0$,\n", "then the Markovh chain is reducible into two recurrent classes, `[1, 2]` and `[3]`:" ] }, { "cell_type": "code", "execution_count": 76, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 0.5 0.5 0.0\n", " 0.5 0.5 0.0\n", " 0.0 0.0 1.0" ] }, "execution_count": 76, "metadata": {}, "output_type": "execute_result" } ], "source": [ "P_epsilon(0)" ] }, { "cell_type": "code", "execution_count": 77, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Array{Array{Int64,1},1}:\n", " [1, 2]\n", " [3] " ] }, "execution_count": 77, "metadata": {}, "output_type": "execute_result" } ], "source": [ "recurrent_classes(MarkovChain(P_epsilon(0)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If $\\varepsilon > 0$ but small, the chain is irreducible,\n", "but transition within each of the subsets `[1, 2]` and `[3]` is much more likely\n", "than that between these sets." ] }, { "cell_type": "code", "execution_count": 78, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 0.499 0.5 0.001\n", " 0.5 0.499 0.001\n", " 0.001 0.001 0.998" ] }, "execution_count": 78, "metadata": {}, "output_type": "execute_result" } ], "source": [ "P_epsilon(0.001)" ] }, { "cell_type": "code", "execution_count": 79, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1-element Array{Array{Int64,1},1}:\n", " [1, 2, 3]" ] }, "execution_count": 79, "metadata": {}, "output_type": "execute_result" } ], "source": [ "recurrent_classes(MarkovChain(P_epsilon(0.001)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Analytically, the unique stationary distribution of the chain with $\\varepsilon > 0$\n", "is `(1/3, 1/3, 1/3)`, independent of the value of $\\varepsilon$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However,\n", "for such matrices with small values of $\\varepsilon > 0$,\n", "general purpose eigenvalue solvers are numerically unstable.\n", "\n", "For example, if we use [`Base.LinAlg.eig`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.eig)\n", "to compute the eigenvector that corresponds\n", "to the dominant (i.e., largest in magnitude) eigenvalue:" ] }, { "cell_type": "code", "execution_count": 80, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "epsilon = 1.0e-11\n", "[0.333333, 0.333333, 0.333333]\n", "epsilon = 1.0e-12\n", "[0.333316, 0.333316, 0.333368]\n", "epsilon = 1.0e-13\n", "[0.333304, 0.333304, 0.333393]\n", "epsilon = 1.0e-14\n", "[0.331651, 0.331651, 0.336699]\n", "epsilon = 1.0e-15\n", "[0.297638, 0.297638, 0.404725]\n", "epsilon = 1.0e-16\n", "[0.5, 0.5, -0.0]\n", "epsilon = 1.0e-17\n", "[0.5, 0.5, -0.0]\n" ] } ], "source": [ "epsilons = [10.0^(-i) for i in 11:17]\n", "for eps in epsilons\n", " println(\"epsilon = $eps\")\n", " \n", " w, v = eig(P_epsilon(eps)')\n", " i = indmax(w)\n", " println(v[:, i]/sum(v[:, i]))\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The same applies to [`Base.LinAlg.eigs`](https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.eigs-Tuple{Any}):" ] }, { "cell_type": "code", "execution_count": 81, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "epsilon = 1.0e-11\n", "[0.333334, 0.333334, 0.333332]\n", "epsilon = 1.0e-12\n", "[0.333333, 0.333333, 0.333333]\n", "epsilon = 1.0e-13\n", "[0.333309, 0.333309, 0.333383]\n", "epsilon = 1.0e-14\n", "[0.333185, 0.333185, 0.33363]\n", "epsilon = 1.0e-15\n", "[0.332257, 0.332257, 0.335485]\n", "epsilon = 1.0e-16\n", "[0.300527, 0.300527, 0.398946]\n", "epsilon = 1.0e-17\n", "[-0.852598, -0.852598, 2.7052]\n" ] } ], "source": [ "epsilons = [10.0^(-i) for i in 11:17]\n", "for eps in epsilons\n", " println(\"epsilon = $eps\")\n", " \n", " w, v = eigs(P_epsilon(eps)', nev=2)\n", " i = indmax(w)\n", " println(v[:, i]/sum(v[:, i]))\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The output becomes farther from the actual stationary distribution `(1/3, 1/3, 1/3)`\n", "as $\\varepsilon$ becomes smaller." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`MarkovChain` in `quantecon` employs\n", "the algorithm called the \"[GTH algorithm](http://quanteconpy.readthedocs.org/en/stable/markov/gth_solve.html)\",\n", "which is a numerically stable variant of Gaussian elimination,\n", "specialized for Markov chains." ] }, { "cell_type": "code", "execution_count": 82, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "epsilon = 1.0e-12\n", "[0.333333, 0.333333, 0.333333]\n", "epsilon = 1.0e-13\n", "[0.333333, 0.333333, 0.333333]\n", "epsilon = 1.0e-14\n", "[0.333333, 0.333333, 0.333333]\n", "epsilon = 1.0e-15\n", "[0.333333, 0.333333, 0.333333]\n", "epsilon = 1.0e-16\n", "[0.333333, 0.333333, 0.333333]\n", "epsilon = 1.0e-17\n", "[0.333333, 0.333333, 0.333333]\n", "epsilon = 1.0e-100\n", "[0.333333, 0.333333, 0.333333]\n" ] } ], "source": [ "epsilons = [10.0^(-i) for i in 12:17]\n", "push!(epsilons, 1e-100)\n", "for eps in epsilons\n", " println(\"epsilon = $eps\")\n", " println(stationary_distributions(MarkovChain(P_epsilon(eps)))[1])\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It succeeds in obtaining the correct stationary distribution for any value of $\\varepsilon$." ] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.6.0", "language": "julia", "name": "julia-0.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.6.0" } }, "nbformat": 4, "nbformat_minor": 1 }