You can use the following websites to run your (R, julia, Rust) code online, without having to install anything:

Slice sampler

This part is taken from Robert and Casella (2004).

At iteration t, simulate:

\begin{eqnarray} u^{(t+1)} &\sim& U_{[0,f(x^{(t)}]} \nonumber \\ x^{(t+1)} &\sim& U_{A^{(t+1)}} \nonumber \end{eqnarray} with \(A^{(t+1)} = \{ x: f(x) \geq u^{(t+1)} \}\)

Example 1

In this example, the slice sampler will be used to simulate from the density: \[ f(x)=\frac{1}{2}\exp( -\sqrt{x} ) \]

At iteration t, simulate: \begin{eqnarray} u^{(t+1)} &\sim& U_{[0,\frac{1}{2}\exp( -\sqrt{x} )]} \nonumber \\ x^{(t+1)} &\sim& U_{[0,(-\log(2u))^2]} \nonumber \end{eqnarray}

Using R

# Sample from 0.5*exp(-sqrt(x))
# using the slice sampler

f <- function(y)exp(-sqrt(y))/2;

nsim <- 50000;

x <- array(0,dim=c(nsim,1));
x[1] <- -log(runif(1));

for (i in 2:nsim)  {
     w <- runif(1,min=0,max=f(x[i-1]));
     x[i] <- runif(1,min=0,max=(-log(2*w))^2);
}

hist(x,main="Slice sample",xlim=c(0,60),ylim=c(0,.25),freq=F,col="green",breaks=250)

par(new=T)

plot(function(y)f(y), 0, 70, xlim=c(0,60),ylim=c(0,.25),xlab="",ylab="",xaxt="n",yaxt="n",lwd=3)

Using julia

using Random, Distributions

# Setting the seed
Random.seed!(123)

function f(y)
    exp( -sqrt(y) ) / 2
end

nsim = 50000

x = rand(Uniform(0.0,1.0),nsim,1)

x[1] = -log.( rand(Uniform(0.0,1.0)) )

for i in 2:nsim
    w = rand(Uniform(0.0,f(x[i-1])))
    x[i] = rand(Uniform(0.0,(-log.(2*w))^2))
end

# Histogram of density
using Plots; pyplot()
Plots.PyPlotBackend()
histogram(x)

Example 2

Now we consider a truncated normal distribution \(\mathcal{N}\)(3,1), which is restricted to the interval [0,1].

In this example, the slice sampler will be used to simulate from the density:

\[f(x) \varpropto f_1(x) = \exp\big( -(x+3)^2/2 \big) \mathbb{1}_{[0,1]}{(x)}\]

At iteration t, simulate: \begin{eqnarray} u^{(t+1)} &\sim& \exp\big( -(x+3)^2/2 \big) U_{[0,1]} \nonumber \\ x^{(t+1)} &\sim& U_{[0,(-3+\sqrt{-2\log(u)})]} \nonumber \end{eqnarray}

# Sample from a N(-3,1) truncated to [0,1]
# using the slice sampler

nsim = 1000

obs <- matrix(0,nsim,2)
obs[1,1] = .25
obs[1,2] = 0.5*exp(-0.5*(obs[1,1]+3)^2)

for (i in 2:nsim){
 obs[i,2] = exp(-0.5*(obs[i-1,1]+3)^2)*runif(1)
 obs[i,1] = runif(1)*min(1,-3+sqrt(-2*log(obs[i,2])))
 }

hist(obs,main="Slice sample",freq=F,col="green", breaks=25)

Monte Carlo Integration

Consider the function, \begin{eqnarray} h(x) = [ \cos(50x) + \sin(20x) ]^2 \nonumber \end{eqnarray}

The integral, \begin{eqnarray} \int^1_0 h(x) dx \nonumber \end{eqnarray} can be calculated by generating \[ U_1, U_2, \ldots, U_n \sim iid \quad U(0,1) \] random variables, and approximate \(\int h(x) dx\) with \[ \sum^n_{i=1} h(U_i)/n \]

Using R

nsim <- 10000;
u <- runif(nsim);
den <- 1:(nsim)

# The function to be integrated
mci.ex <- function(x){(cos(50*x)+sin(20*x))^2}
par(mfrow=c(1,3))
plot(function(x)mci.ex(x), xlim=c(0,1),ylim=c(0,4),xlab="Function",ylab="")

# The Monte Carlo sum
## Generated Values of Function
hint <- mci.ex(u)
## Mean
hplot <- cumsum(hint)/den
## Standard Errors
stdh <- sqrt( cumsum(hint^2)/den - (hplot)^2)
upper <- hplot+stdh/sqrt(den)
lower <- hplot-stdh/sqrt(den)

Simulated Annealing

Given a temperature parameter \(T>0\), a sample \(\theta^T_1, \theta^T_2, \ldots\) is generated from the distribution \[ \pi(\theta) \varpropto \exp( h(\theta)/T ) \]

Example 3

We are searching for the maximum of the following function \[ h(x) = [ \cos(50x) + \sin(20x) ]^2 \] Note: Variable

At iteration t the algorithm is at \(( x^{(t)},h^{(t)} )\):

  1. Simulate \(u \sim U(a_t, b_t)\) where \(a_t = \max( x^{(t)}-r,0 )\) and \(b_t = \min( x^{(t)}+r,1 )\)

  2. Accept \(x^{(t+1)}=u\) with probability \[ \rho^{(t)} = \min \big\{ \exp \big( \frac{h(u)-h(x^{(t)})}{T_t} \big), 1 \big\} \] take \(x^{(t+1)}=x^{(t)}\) otherwise.

  3. Update \(T_t\) to \(T_{t+1}\)

par(mfrow=c(1,2))

# The function to be optimized
mci <- function(x){(cos(50*x)+sin(20*x))^2}

# The Monte Carlo maximum
nsim <- 2500
u <- runif(nsim)

# Simulated annealing
xval <- array(0,c(nsim,1));
r <- .5
for(i in 2:nsim){
    test <- runif(1, min=max(xval[i-1]-r,0),max=min(xval[i-1]+r,1));
    delta <- mci(test) - mci(xval[i-1]);
    rho <- min(exp(delta*log(i)/1),1);
    xval[i] <- test*(u[i]<rho)+xval[i-1]*(u[i]>rho)
}
mci(xval[nsim])

# Plot the trajectory of the optimization path
plot(function(x)mci(x), xlim=c(0,1),ylim=c(0,4),xlab="Function",ylab="")
plot(xval,mci(xval),type="l",lwd=2)

Metropolis-Hastings Algorithm

Let \(q(\theta, \vartheta)\) be a proposal density or a candidate-generating density (Chib and Greenberg, 1995) such that \[ \int q(\theta, \vartheta) d\vartheta = 1 \] Also let \(U(O, 1)\) denote the uniform distribution over \((0, 1)\). Then, a general version of the Metropolis- Hastings algorithm for sampling from the posterior distribution \(\pi(\theta,D)\) can be described as follows:

  1. Choose an arbitrary starting point \(\theta_0\) and set \(i=0\).

  2. Generate a candidate point \(\theta^*\) from \(q(\theta_i,\cdot)\) and \(u\) from \(U(0,1)\).

  3. Set \(\theta_{i+1}=\theta^*\) if \(u \leq a(\theta_i, \theta^*)\) and \(\theta_{i+1}=\theta_i\) otherwise, where the acceptance probability is given by \[ a(\theta,\vartheta) = \min\big\{ \frac{ \pi(\vartheta|D)q(\vartheta, \theta) }{ \pi(\theta|D)q(\theta, \vartheta) }, 1 \big\} \]

  4. Set \(i=i+1\), and go to Step 2.

Example 4

Let us consider the following simple linear model: \[ y_t = \alpha x_t + \beta + \epsilon_t \] where \(\epsilon \sim N(0, \sigma^2)\).

Source: R code taken from A simple Metropolis-Hastings MCMC in R

### Creating test data
trueA <- 5
trueB <- 0
trueSd <- 10
sampleSize <- 31

# create independent x-values
x <- (-(sampleSize-1)/2):((sampleSize-1)/2)
# create dependent values according to ax + b + N(0,sd)
y <-  trueA * x + trueB + rnorm(n=sampleSize,mean=0,sd=trueSd)

plot(x,y, main="Test Data")


### Deriving the likelihood
likelihood <- function(param){
    a = param[1]
    b = param[2]
    sd = param[3]

    pred = a*x + b
    singlelikelihoods = dnorm(y, mean = pred, sd = sd, log = T)
    sumll = sum(singlelikelihoods)
    return(sumll)
}

# Example: plot the likelihood profile of the slope a
slopevalues <- function(x){return(likelihood(c(x, trueB, trueSd)))}
slopelikelihoods <- lapply(seq(3, 7, by=.05), slopevalues )
plot (seq(3, 7, by=.05), slopelikelihoods , type="l", xlab = "values of slope parameter a", ylab = "Log likelihood")


### Defining the prior
# Prior distribution
prior <- function(param){
    a = param[1]
    b = param[2]
    sd = param[3]
    aprior = dunif(a, min=0, max=10, log = T)
    bprior = dnorm(b, sd = 5, log = T)
    sdprior = dunif(sd, min=0, max=30, log = T)
    return(aprior+bprior+sdprior)
}

### Defining the posterior
posterior <- function(param){
   return (likelihood(param) + prior(param))
}



######## Metropolis algorithm ################

proposalfunction <- function(param){
    return(rnorm(3,mean = param, sd= c(0.1,0.5,0.3)))
}

run_metropolis_MCMC <- function(startvalue, iterations){
    chain = array(dim = c(iterations+1,3))
    chain[1,] = startvalue
    for (i in 1:iterations){
        proposal = proposalfunction(chain[i,])

        probab = exp(posterior(proposal) - posterior(chain[i,]))
        if (runif(1) < probab){
            chain[i+1,] = proposal
        }else{
            chain[i+1,] = chain[i,]
        }
    }
    return(chain)
}

startvalue = c(4,0,10)
chain = run_metropolis_MCMC(startvalue, 10000)

burnIn = 5000
acceptance = 1-mean(duplicated(chain[-(1:burnIn),]))


### Summary: #######################

par(mfrow = c(2,3))
hist(chain[-(1:burnIn),1],nclass=30, , main="Posterior of a", xlab="True value = red line" )
abline(v = mean(chain[-(1:burnIn),1]))
abline(v = trueA, col="red" )
hist(chain[-(1:burnIn),2],nclass=30, main="Posterior of b", xlab="True value = red line")
abline(v = mean(chain[-(1:burnIn),2]))
abline(v = trueB, col="red" )
hist(chain[-(1:burnIn),3],nclass=30, main="Posterior of sd", xlab="True value = red line")
abline(v = mean(chain[-(1:burnIn),3]) )
abline(v = trueSd, col="red" )
plot(chain[-(1:burnIn),1], type = "l", xlab="True value = red line" , main = "Chain values of a", )
abline(h = trueA, col="red" )
plot(chain[-(1:burnIn),2], type = "l", xlab="True value = red line" , main = "Chain values of b", )
abline(h = trueB, col="red" )
plot(chain[-(1:burnIn),3], type = "l", xlab="True value = red line" , main = "Chain values of sd", )
abline(h = trueSd, col="red" )

# for comparison:
summary(lm(y~x))

The following table gives the equivalence of the variables used in the pseudo-code to the variables used in the R code:

Theory Code

\(\theta\)

chain[,]

\(u\)

runif(1)

\(a(\theta_i, \theta^*)\)

probab

\(\min \big\{ \frac{\pi(\vartheta|D)q(\vartheta,\theta)}{\pi(\theta|D)q(\theta,\vartheta)} \big\}\)

exp( posterior(proposal) - posterior(chain[i,]) )

\(\theta_0\)

c(4,0,10)

\(\theta^*\)

proposal

Dynamic programming and value function iteration

Consider the following dynamic optimization problem: \[ \max_{ \{ x_{t+1} \}^{\infty}_{t=0} } \Sigma^{\infty}_{t=0}\beta^t F(x_t, x_{t+1}) \\ \mbox{s.t.} \quad x_{t+1} \in \Gamma(x_t) \]

The functional equation of the previous problem is given by: \[ v(x) = \max_{y\in \Gamma(x)} \big[ F(x,y) + \beta v(y) \big] \]

Introducing the optimal policy function \(g\), we can rewrite the previous expression: \[ v(x) = F\big[ x, g(x) + \beta v[g(x)] \big] \]

The first-order condition and the envelope condition is defined as: \begin{eqnarray} 0 &=& F_y[ x, g(x) ] + \beta v'[g(x)] \nonumber \\ v'(x) &=& F_x[ x, g(x) ] \nonumber \end{eqnarray}

The Euler equation and the transversality condition is given by: \begin{eqnarray} 0 &=& F_y(x_t, x_{t+1}) + \beta F_x (x_{t+1}, x_{t+2}), \quad t=0,1,2, \ldots \nonumber \\ 0 &=& \lim_{t\rightarrow \infty} \beta^t F_x(x_t, x_{t+1}) \cdot x_t \nonumber \end{eqnarray}

Example 5

The Model

Let us consider the following functional equation: \[ v(k) = \max_{k'\geq 0} \big\{ U\big( f(k)-k' \big) + \beta v(k') \big\} \] with \(k'+c\leq f(k)\), \(c\geq 0\) and \(y=f(k)\)

The model is specified as follows: \[ f(k) = k^{\alpha} \\ U(c) = \log(c) \] and it is parameterized as follows \begin{eqnarray} \alpha &=& 0.3 \nonumber \\ \beta &=& 0.9 \nonumber \end{eqnarray}

First Order and Envelope Conditions

\begin{eqnarray} v(k) &=& \max_{c,k'} \big\{ U(c) + \beta v(k') \big\} \nonumber \\ && \mbox{s.t.} \quad f(k) - k' - c \geq 0 \quad \quad [\lambda] \nonumber \end{eqnarray}

\begin{eqnarray} \frac{\partial v(k)}{\partial k'} = 0 \iff \beta \frac{\partial v(k')}{\partial k'} - \lambda = 0 \nonumber \end{eqnarray}

\begin{eqnarray} \boxed{\lambda = \beta \frac{\partial v(k')}{\partial k'} } \nonumber \end{eqnarray}

\begin{eqnarray} \frac{\partial v(k)}{\partial k} &=& \beta \frac{\partial v(k')}{\partial k'} + \lambda \big[ \frac{\partial f(k)}{\partial k} - \frac{k'}{k} \big] \nonumber \\ &=& \underbrace{\big[ \beta \frac{\partial v(k')}{\partial k'} - \lambda \big]}_{=0} \frac{k'}{\partial k} + \lambda \frac{\partial f(k)}{\partial k} \nonumber \end{eqnarray} which yields the envelope condition \[ \boxed{ \frac{\partial v(k')}{\partial k'} = \lambda' \frac{\partial f(k')}{\partial k'} } \]

By combining the two boxed expressions, we get the first-order condition \[ \boxed{ \lambda = \lambda' \beta \frac{\partial f(k')}{\partial k'} } \] and \[ \frac{\partial v(k)}{\partial c} = 0 \iff \frac{\partial U(c)}{\partial c} - \lambda = 0 \]

\[ \boxed{ \lambda = \frac{\partial U(c)}{\partial c} } \]

Steady-state

By combining the first-order condition with the following expressions \begin{eqnarray} \frac{\partial U(c)}{\partial c} &=& c^{-1} \nonumber \\ \frac{\partial f(k)}{\partial k} &=& \alpha k^{\alpha-1} \nonumber \end{eqnarray} we can find the steady-state value of \(k\) \[ k^* = \Bigg( \frac{1}{\alpha\beta} \Bigg)^{\frac{1}{\alpha-1}} \]

Numerical Solution

The numerical solution algorithm (value function iteration) is given by (Heer and Maussner, 2008):

  1. Choose a grid \(\mathcal{K} = [ K_1, K_2, \ldots, K_n\)] with \(K_i < K_j\) and \(i < j = 1\),\(2\),\(\ldots\), \(n\)

  2. Initialize the value function: \(\forall i = 1, \ldots, n\) set \[ v_i^0 = \frac{u(f(K^*)-K^*)}{1-\beta} \] where \(K^*\) denotes the stationary solution to the initial problem

  3. Compute a new value function and the associated policy function, \(\underline{v}^1\) and \(\underline{h}^1\) respectively: Put \(j^*_0 \equiv 1\). For \(i=1,2,\ldots , n\) and \(j^*_{i-1}\) find the index \(j^*_i\) that maximizes \[ u\big( f(K_i)-K_j \big) + \beta v^0_j \] in the set of indices \(\{ j^*_{i-1}, j^*_{i-1}+1, \ldots, n \}\). Set \(h^1_i=j^*_i\) and \(v^1_i=u\big( f(K_i)-K_{j^*_i} \big) + \beta v^0_{j^*_i}\)

  4. Check for convergence: If \(\parallel \underline{v}^0-\underline{v}^1 \parallel_{\infty} < \epsilon (1-\beta)\), \(\epsilon \in \mathbb{R}_{++}\) stop, else replace \(\underline{v}^0\) with \(\underline{v}^1\) and \(\underline{h}^0\) with \(\underline{h}^1\) and return to step 3.

Using Octave

% Deterministic Growth model
% Original code by Joao Ejarque and Russell Cooper
% Tarik OCAKTAN - Paris - May 2009

% State space dimension = 1

% Notation:
%  pName : structural parameters
%  iName : total number of variables related to the algorithm
%  idxName : indexation used in loops, i.e. idxName = 1,...,iName
%  dName : integer of a specific value of a variable, i.e. dNamess=steady state
%          of
%							   variable 'Name'
%  vName : vector containing values of variable 'Name'
%  mName : matrix containing values of a variable 'Name'

clear;clc;
close all;

% Structural parameters
pAlpha	= .3;
pBeta	= .9;

% Algorithm parameters
 iIter = 2000; % maximum number of iterations
 iGridPoints = 300;
 iToler	= 10e-6;  % convergence criterion
% Steady state value
 dKss = ( 1 / (pAlpha*pBeta) )^( 1 / (pAlpha-1) );
% Bounds of state space
 dKmin	= dKss * 0.9;
 dKmax	= dKss * 1.1;
% Construct a grid
 iStep	= (dKmax-dKmin) / (iGridPoints-1);
 vK	= dKmin:iStep:dKmax;
 iGridPoints = size(vK,2); % overwriting iGridPoints in order to take into
			   % account additional point

vY = vK.^pAlpha;
vIota = ones(iGridPoints,1);

% Budget constraint
mC = vIota*vY - vK'*vIota';

% Eliminate negative or zero consumption values
mCP	= mC > 0;  % Matrix of ones where C>0 and zeros elsewhere
mPC	= 1 - mCP; % Matrix with ones where C<=0. this is a penalty matrix
z	= 10e-5;
mC	= mC.*mCP + z*mPC;

% Initialization
 % Initializing the utility function mU
 mU = log(mC);

 % Initializing the value function vVNew
 mCF = vIota*vY ;
 mVNew = log(mCF);

 % Initializing the value function vVOld
 vVOld = max([ mU + pBeta*mVNew ]);

%%% Main Iteration - BEGIN %%%
for idxIter=1:iIter

 idxIter

 vVNew = max([ mU + pBeta*(vIota*vVOld)' ]);

 iDiff = ( abs(vVOld-vVNew) ) / vVOld;
 if abs(iDiff) <= iToler
  vVOld = vVNew;
  break
 else
  vVOld = vVNew;
 end

end
%%% Main Iteration - END %%%



%break

% Construct policy function
[maxR,idxR] = max([ mU + pBeta*(vIota*vVOld)' ]);

vKP  = vK(idxR);
vInv = vKP;

figure(1)
plot(vK,vKP)
figure(2)
plot(vK,vY)
disp(idxIter)

Alternative version TBD

% Deterministic Growth model
% Original code by Eric Sims
% Tarik OCAKTAN - Luxembourg - November 2020

% State space dimension = 1

% Notation:
%  pName : structural parameters
%  iName : total number of variables related to the algorithm
%  idxName : indexation used in loops, i.e. idxName = 1,...,iName
%  dName : integer of a specific value of a variable, i.e. dNamess=steady state
%          of
%							   variable 'Name'
%  vName : vector containing values of variable 'Name'
%  mName : matrix containing values of a variable 'Name'

clear;clc;
close all;

global v0 pBeta pDelta pAlpha vK k0 s

% Structural parameters
pAlpha	=  .33;
pBeta	=  .95;
pDelta  = 0.1;
s       = 2;

% Value function
function val=valfun(k)

    % This program gets the value function for a neoclassical growth model with
    % no uncertainty and CRRA utility
    global v0 pBeta pDelta pAlpha vK k0 s
    g = interp1(vK,v0,k,'linear'); % smooths out previous value function
    c = k0^pAlpha-k + (1-pDelta)*k0; % consumption
    if c<=0
        val =-888888888888888888-800*abs(c); % keeps it from going negative
    else
        val=(1/(1-s))*(c^(1-s)-1) + pBeta*g;
    end
    val =-val; % make it negative since we're maximizing and code is to minimize.

endfunction


% Algorithm parameters
 iIter = 1000; % maximum number of iterations
 iGridPoints = 4;
 iToler	= 0.01;  % convergence criterion

% Steady state value
 dKss = (pAlpha/(1/pBeta-(1-pDelta)))^(1/(1-pAlpha));
 dCss = dKss^(pAlpha)-pDelta*dKss;
 dIss = pDelta*dKss;
 dYss = dKss^pAlpha;

% Bounds of state space
 dKmin	= dKss * 0.25;
 dKmax	= dKss * 1.75;
% Construct a grid
 iStep	= (dKmax-dKmin) / (iGridPoints-1);
 vK	= [dKmin:iStep:dKmax]';
 iGridPoints = size(vK,1); % overwriting iGridPoints in order to take into
			   % account additional point

 polfun = zeros(iGridPoints,1);

 v0 = zeros(iGridPoints,1);

 iDif = 10;
 idxIter = 0;



while iDif>iToler && idxIter<iIter

    for i = 1:iGridPoints

        k0 = vK(i,1);
        k1 = fminbnd(@valfun,dKmin,dKmax);
        v1(i,1) = -valfun(k1);
        k11(i,1) = k1;

    end

iDif = norm(v1-v0);
v0 = v1;
idxIter = idxIter+1;

end

break

for i = 1:iGridPoints
    con(i,1) = vK(i,1)^(pAlpha)-k11(i,1) + (1-pDelta)*vK(i,1);
    polfun(i,1) = vK(i,1)^(pAlpha)-k11(i,1) + (1-pDelta)*vK(i,1);
end

Alternative version This julia code contains an error and will be fixed soon …​

Example 6

Let us consider the following two-dimensional deterministic model: \begin{eqnarray} v(k,n) = \max_{ n',k',c } \Bigg\{ U(c) - G(n) &+& \beta v(k',n') \Bigg\} \nonumber \\ \mbox{ s.t. } \quad f(k,n)-k'-c-(1-n) &\geq& 0 \quad [\lambda] \nonumber \nonumber \\ \quad \quad \quad (1-s)n + p(1-n) - n' &\geq& 0 \quad [\mu] \nonumber \end{eqnarray}

The first-order conditions describing the dynamics of this system are given by:

\begin{eqnarray} k^{\theta} n^{1-\theta} - k' - c - 1 + n &=& 0 \nonumber \\ (1-s)n + p(1-n) - n' &=& 0 \nonumber \end{eqnarray}

\[ \beta (c')^{-1} \theta k^{ \theta - 1 } n^{ 1 - \theta } = c^{ -1 } \]

\[ \beta \Bigg[ -(n')^{-1} + (c')^{-1} \Bigg( 1+(1-\theta)k^{\theta}n^{-\theta} \Bigg) + \mu'(1-s-p) \Bigg] = \mu \nonumber \]

Using the previous equations, we can compute the steady state values for the four variables, which is given by: \[ k^* = \Bigg[ \frac{1}{\beta\theta\big( \frac{p}{s+p} \big)^{1-\theta}} \Bigg]^{\frac{1}{\theta-1}} \\ n^* = \frac{p}{s+p} \\ c^* = (k^*)^{\theta} (n^*)^{1-\theta} - k^* - 1 + n^* \\ \mu^* = \frac{ \beta\big[ -(n^*)^{-1} + (c^*)^{-1} \big( 1+(1-\theta)(k^*)^{\theta}(n^*)^{-\theta} \big) \big] }{1-\beta(1-s-p)} \]

The model is parameterized as follows:

Parameter Value

\(\beta\)

0.95

\(\theta\)

0.2

\(s\)

0.3

\(p\)

0.75

This model is a modified version of Merz, 1995 and Andolfatto, 1996.

Using Octave

clear;clc;
close all;

% Structural parameters
Parameters.pBeta = .95;
Parameters.pTheta = .2;
Parameters.pS = .3;
Parameters.pP = .75;

% Algorithm parameters
 iIter = 2000; % maximum number of iterations
 Grid.nK = 3;
 Grid.nN = 2;
 iToler	= 10e-6;  % convergence criterion

% Steady state value
 dKss = ( 1 / (Parameters.pBeta*Parameters.pTheta*(Parameters.pP/(Parameters.pS+Parameters.pP))^(1-Parameters.pTheta)) )^( 1 / (Parameters.pTheta-1) );
 dNss = Parameters.pP/(Parameters.pS+Parameters.pP);
 dCss = (dKss)^(Parameters.pTheta)*(dNss)^(1-Parameters.pTheta)-dKss-1+dNss;
 dMu = ( Parameters.pBeta*( -dNss^(-1)+dCss^(-1)*( 1+(1-Parameters.pTheta)*dKss^Parameters.pTheta * dNss^(-Parameters.pTheta) ) ) )/( 1-Parameters.pBeta*(1-Parameters.pS-Parameters.pP) );

% Bounds of state space
 dKmin	= dKss * 0.9;
 dKmax	= dKss * 1.1;
 dNmin = dNss * 0.9;
 dNmax = dNss * 1.1;

% Create a grid for K and N
Grid.K = linspace(dKmin, dKmax,Grid.nK)';
Grid.N = linspace(dNmin, dNmax,Grid.nN)';

% Storage space for consumption
mC = ones( (Grid.nN-1)*Grid.nK , Grid.nK);
mY = ones( Grid.nN, Grid.nK );

for iN=1:Grid.nN
    for iKp=1:Grid.nK
        for iK=1:Grid.nK
            mY( iN, iK ) = Grid.K(iK)^Parameters.pTheta*Grid.N(iN)^(1-Parameters.pTheta);
            mC( (iN-1)*Grid.nK+iKp, iK ) = Grid.K(iK)^Parameters.pTheta*Grid.N(iN)^(1-Parameters.pTheta) - Grid.K(iKp) - 1 + Grid.N(iN);
        end
    end
end

% Eliminate negative or zero consumption
cp = mC > 0;
pc = 1 - cp;
z = 10e-5;
mC = mC.*cp + z*pc;

% Initialize the utility function
mU = log( mC );


% Initializing the value function vVNew
mCF = ones( Grid.nN*Grid.nK, Grid.nN )*mY ;
mVNew = log(mCF);


% Initializing the value function vVOld
vVOld = max([ mU + Parameters.pBeta*mVNew ]);


%%% Main Iteration - BEGIN %%%
for idxIter=1:iIter

 idxIter

 vVNew = max([ mU + Parameters.pBeta*( ones( Grid.nN*Grid.nK, 1 )*vVOld ) ]);

 iDiff =  ( abs(vVOld-vVNew) ) / vVOld ;


 if abs(iDiff) <= iToler
  vVOld = vVNew;
  break
 else
  vVOld = vVNew;
 end

end
%%% Main Iteration - END %%%


% Construct policy function
[maxR,idxR] = max([ mU + Parameters.pBeta*( ones( Grid.nN*Grid.nK, 1 )*vVOld ) ]');

maxR
idxR

break

% INCORRECT - TBD

vKP  = Grid.K(idxR);
vInv = vKP;

figure(1)
plot(Grid.K,vKP(1:Grid.nK))
print -dpng some_name.png;


figure(1)
plot(Grid.K,vKP(1:Grid.nK))

break
figure(2)
plot(Grid.K,mY(1,:))
disp(idxIter)

Example 7

Let us consider the following two-dimensional deterministic model: \begin{eqnarray} w(k,n) = \max_{ c, v, s, k', n' } \Bigg\{ u(c) - g(n) && \beta v(k',n') \Bigg\} \nonumber \\ \mbox{ s.t. } \quad f(k,n) - k' + (1\delta) -c(s)\cdot (1-n) - \gamma v - c &\geq& 0 \quad [\lambda] \nonumber \nonumber \\ \quad \quad \quad (1-\rho_x)n + m(n,v) - n' &\geq& 0 \quad [\mu] \nonumber \end{eqnarray}

where \(u(c,n)=u(c)-g(n)\) with \(u(c)=\log(c)\) and \(g(n)= \theta \frac{n^{1-\gamma_n}}{1-\gamma_n}\) ; \(m(n.v)=\mu_m v^{1-\alpha_m}[s(1-n)^{\alpha_m}\)] ; \(c(s)=c_0 s^{\eta}\) ; \(f(k,n)=k^{\alpha}n^{1-\alpha}\) .

The model is parameterized as follows:

Parameter Value

\(\alpha\)

0.36

\(\alpha_m\)

0.40

\(\theta\)

1

\(\gamma_n\)

\(-\frac{1}{1.25}\)

\(c_0\)

0.005

\(\delta\)

0.022

\(\rho_x\)

0.07

\(\mu_m\)

1

\(\gamma\)

0.05

\(\beta\)

0,99

The steady state values for this model are:

Parameter Value

\(c\)

2.721

\(v\)

0.039

\(s\)

1.917

\(k\)

40.553

\(n\)

0.928

The following code is an implementation of the perturbation method for solving this model.

using Dynare

var c k n v i s f f_k f_n m m_v m_s m_n cs u_c g_n cs_s ;

varexo z;

parameters  	pAlpha pBeta pDelta pEta pGamma_n pAlpha_m
		pRho_x pMu pA pC0 pRho pSigmaeps pMuz;

pAlpha	= 0.36;
pBeta 	= 0.99;
pDelta 	= 0.022;
pEta 	= 1;
pGamma_n 	= -1/1.25;
pAlpha_m	= 0.40;
pRho_x 	= 0.07;
pMu 	= 0.004;
pA 	= 0.05;
pC0 	= 0.005;
pRho 	= 0.95;
pSigmaeps = 0.007;

pMuz	= 1;

model;
// Production function
f = z*k(-1)^pAlpha*n(-1)^(1-pAlpha);
// Marginal productivity of capital
f_k = pAlpha*z*k(-1)^(pAlpha-1)*n(-1)^(1-pAlpha);
// Marginal productivity of labor
f_n = (1-pAlpha)*z*k(-1)^pAlpha*n(-1)^(-pAlpha);
// Matching function
m = v^(1-pAlpha_m)*( s*( 1-n(-1) ) )^pAlpha_m;
// Matching function - derivative w.r.t. vacancy
m_v = (1-pAlpha_m)*v^(-pAlpha_m)*( s*( 1-n(-1) ) )^pAlpha_m;
// Matching function - derivative w.r.t. search effort
m_s = pAlpha_m*( 1-n(-1) )*v^(1-pAlpha_m)*( s*(1-n(-1)) )^(pAlpha_m-1);
// Matching function - derivative w.r.t. employment
m_n = -pAlpha_m*s*v^(1-pAlpha_m)*( s*( 1-n(-1) ) )^(pAlpha_m-1);
// Marginal utility for consumption
u_c = 1/c;
// Marginal utility for labor
g_n = n(-1)^(-1/pGamma_n);
// Cost for search effort
cs = pC0*s^pEta;
// Cost for search effort - derivative w.r.t. search effort
cs_s = pEta*pC0*s^(pEta-1);
// P1
u_c = pBeta*( u_c(+1)*( f_k(+1) + (1-pDelta) ) );
// P2
( pA*u_c ) / m_v = pBeta*
	(
	 u_c(+1)*( f_n(+1)+cs(+1) )
	 - g_n(+1)
	 + ( ( pA*u_c(+1) ) / m_v(+1) ) * ( (1-pRho_x)+m_n(+1) )
	);
// P3
( u_c*cs*( 1-n(-1) ) ) / m_s = pBeta*
	(
	 u_c(+1) * ( f_n(+1) + cs(+1) )
	 - g_n(+1)
	 + ( ( u_c(+1)*cs(+1)*( 1-n ) ) / m_s(+1) ) * ( (1-pRho_x)+m_n(+1) )
	);
// Budget constraint
c = f - i - cs*( 1-n(-1) ) - pA*v;
// Investment
k = (1-pDelta)*k(-1) + i;
// Employment dynamics
n = (1-pRho_x)*n(-1) + m;
end;

initval;
c    = 2.67327;
k    = 39.8371;
n    = 0.912013;
z    = 1;
v    = 0.0366684;
i    = 0.876417;
s    = 1.66683;
f    = 3.55226;
f_k  = 0.032101;
f_n  = 2.49278;
m    = 0.0638409;
m_v  = 1.04462;
m_s  = 0.0153203;
m_n  = -0.290228;
cs   = 0.00833414;
u_c  = 0.374073;
g_n  = 0.928968;
cs_s =  0.005;
end;

steady;

endval;
 z = 1;
end;

check;

// Declare a positive technological shock in period 1
//shocks;
//var z; stderr 0.001;
//end;


//stoch_simul(order=1);
simul(periods=300);

The following code gives an implementation of the value function iteration. For this we need to construct a grid for the state space of the model, which is givven in the following diagram:

VFIGRID

using Octave

%pAlpha, pAlpha_m, pTheta, pGamma_n, pC0, pEta, pDelta, pRho_x, pMu_m, pGamma, pBeta

% Structural parameters
global pAlpha = 0.36;
global pAlpha_m = 0.40;
global pTheta = 1;
global pGamma_n = -1/1.25;
global pC0 = 0.005;
global pEta = 1;
global pDelta = 0.022;
global pRho_x = 0.07;
global pMu_m = 1;
global pGamma = 0.05;
global pBeta = 0.99;

% Computing the steady state

PARAMS = [pAlpha, pAlpha_m, pTheta, pGamma_n, pC0, pEta, pDelta, pRho_x, pMu_m, pGamma, pBeta];

% 1  2  3  4  5    6     7
% c, v, s, k, n, lambda, mu
function y = f(x, PARAMS)

    pAlpha = PARAMS(1);
    pAlpha_m = PARAMS(2);
    pTheta = PARAMS(3);
    pGamma_n = PARAMS(4);
    pC0 = PARAMS(5);
    pEta = PARAMS(6);
    pDelta = PARAMS(7);
    pRho_x = PARAMS(8);
    pMu_m = PARAMS(9);
    pGamma = PARAMS(10);
    pBeta = PARAMS(11);


     y(1) = x(4)^pAlpha*x(5)^(1-pAlpha) -pDelta*x(4) -pC0*x(3)^pEta*(1-x(5))-pGamma*x(2)-x(1);
     y(2) = (1-pRho_x)*x(5) + pMu_m*x(2)^(1-pAlpha_m)*( x(3)*(1-x(5)) )^pAlpha_m - x(5);
     y(3) = pBeta*( pAlpha*x(4)^(pAlpha-1)*x(5)^(1-pAlpha) + (1-pDelta) ) - 1;
     y(4) = pBeta*( -pTheta*x(5)^(-pGamma_n) + x(6)*( (1-pAlpha)*x(4)^pAlpha*x(5)^(-pAlpha)+pC0*x(3)^pEta ) + x(7)*( (1-pRho_x)-pMu_m*pAlpha_m*x(2)^(1-pAlpha_m)*x(3)^pAlpha_m*(1-x(5))^(pAlpha_m-1) ) ) - x(7);
     y(5) = 1/x(1) - x(6);
     y(6) = x(7)*( pMu_m*(1-pAlpha_m)*x(2)^(-pAlpha_m)*( x(3)*(1-x(5)) )^(pAlpha_m) )-pGamma*x(6);
     y(7) = x(7)*pMu_m*x(2)^(1-pAlpha_m)*pAlpha_m*x(3)^(pAlpha_m-1)*(1-x(5))^pAlpha_m-x(6)*pC0*pEta*x(3)^(pEta-1)*(1-x(5));

endfunction

%[x, info] = fsolve ("f", [2.5; 0.03; 1.9; 30; 0.8; 1; 0.1; 0.02; 1])
f = @(y) f(y,PARAMS); % function of dummy variable y
[x,fval]=fsolve(f,[2.5; 0.03; 1.9; 30; 0.8; 1; 0.1; 0.02; 1]);



% The steady state values
dCss = x(1);
dVss = x(2);
dSss = x(3);
dKss = x(4);
dNss = x(5);
dLambda_ss = x(6);
dMu_ss = x(7);


% Algorithm parameters
iIter = 2000; % maximum number of iterations
Grid.nK = 30;
Grid.nN = 20;
iToler	= 10e-6;  % convergence criterion

% Bounds of state space
dKmin	= dKss * 0.99;
dKmax	= dKss * 1.01;
dNmin = dNss * 0.99;
dNmax = dNss * 1.01;

% Create a grid for K and N
Grid.K = linspace(dKmin, dKmax,Grid.nK)';
Grid.N = linspace(dNmin, dNmax,Grid.nN)';

% Storage space for consumption
mC = ones( (Grid.nN-1)*Grid.nK , Grid.nK);
mY = ones( Grid.nN, Grid.nK );

for iNp=1:Grid.nN
    for iN=1:Grid.nN
        for iKp=1:Grid.nK
            for iK=1:Grid.nK
                mY( iN, iK ) = Grid.K(iK)^pAlpha*Grid.N(iN)^(1-pAlpha);
                % For given n and n', solve for v and s:
                % 1 2
                % v s
                REF = [Grid.N(iNp), Grid.N(iN), PARAMS];
                function y = functionVS(x, REF)
                  %
                  np = REF(1);
                  n = REF(2);
                  pAlpha = REF(3);
                  pAlpha_m = REF(4);
                  pTheta = REF(5);
                  pGamma_n = REF(6);
                  pC0 = REF(7);
                  pEta = REF(8);
                  pDelta = REF(9);
                  pRho_x = REF(10);
                  pMu_m = REF(11);
                  pGamma = REF(12);
                  pBeta = REF(13);
                  %
                  y(1) = ( (np-(1-pRho_x)*n)*pMu_m^(-1)*( x(2)*(1-n) )^(-pAlpha_m) )^(1/(1-pAlpha_m)) ...
                            - x(1);
                  y(2) =  pC0*pEta*pGamma^(-1)*pAlpha_m^(-1)*(1-pAlpha_m)*(1-n)*x(1)^(-1) - x(2)^(-pEta);
                endfunction
                %[M, info] = fsolve ("functionVS", [0.03;1.9])
                f = @(y) functionVS(y,REF); % function of dummy variable y
                [M,fval]=fsolve(f,[0.03;1.9]);
                dV = M(1);
                dS = M(2);
                % Compute the consumption
                mC( (iN-1)*Grid.nK+iKp, iK ) = Grid.K(iK)^pAlpha*Grid.N(iN)^(1-pAlpha) ...
                                                - Grid.K(iKp) ...
                                                + (1-pDelta)*Grid.K(iK) ...
                                                - pC0*dS^pEta*(1-Grid.N(iN)) ...
                                                - pGamma*dV;
            end
        end
    end
end


% Eliminate negative or zero consumption
cp = mC > 0;
pc = 1 - cp;
z = 10e-5;
mC = mC.*cp + z*pc;


% Initialize the utility function
mU = log( mC );


% Initializing the value function vVNew
mCF = ones( Grid.nN*Grid.nK, Grid.nN )*mY ;
mVNew = log(mCF);


% Initializing the value function vVOld
vVOld = max([ mU + pBeta*mVNew ]);


%%% Main Iteration - BEGIN %%%
for idxIter=1:iIter

 %idxIter

 vVNew = max([ mU + pBeta*( ones( Grid.nN*Grid.nK, 1 )*vVOld ) ]);

 iDiff =  ( abs(vVOld-vVNew) ) / vVOld ;


 if abs(iDiff) <= iToler
  vVOld = vVNew;
  break
 else
  vVOld = vVNew;
 end

end
%%% Main Iteration - END %%%


% Construct policy function
[maxR,idxR] = max([ mU + pBeta*( ones( Grid.nN*Grid.nK, 1 )*vVOld ) ]');

maxR
idxR


disp("done!")

This is julia implementation of a 2D value iteration algorithm (incomplete graph): Julia code for value function iteration of a 2D problem.

Functional approximation: Piecewise linear interpolation

If \((x_i, y_i)\) is the data, then the piecewise linear interpolant \(\hat{f}(x)\) is given by the formula: \begin{eqnarray} \hat{f}(x) = y_i + \frac{x-x_i}{x_{i+1}-x_i} (y_{i+1}-y_i) \quad \quad x\in[x_i,x_{i+1}] \nonumber \end{eqnarray} The julia code for this interpolation is given by:

function interpolate1D(vX, vY, x)
    imin = sum(vX.<x);
    imax = sum(vX.<x)+1;
    y = vY[imin]+( x-vX[imin] )/( vX[imax]-vX[imin] )*(vY[imax]-vY[imin]);
    return y
end