Permalink
Fetching contributors…
Cannot retrieve contributors at this time
108 lines (73 sloc) 3.19 KB
## 0. Housekeeping
rm(list=ls())
ptm <- proc.time()
## 1. Calibration
aalpha = 1/3; # Elasticity of output w.r.t. capital
bbeta = 0.95; # Discount factor
# Productivity values
vProductivity <- c(0.9792, 0.9896, 1.0000, 1.0106, 1.0212);
# Transition matrix
mTransition <- c(0.9727, 0.0273, 0.0000, 0.0000, 0.0000,
0.0041, 0.9806, 0.0153, 0.0000, 0.0000,
0.0000, 0.0082, 0.9837, 0.0082, 0.0000,
0.0000, 0.0000, 0.0153, 0.9806, 0.0041,
0.0000, 0.0000, 0.0000, 0.0273, 0.9727);
mTransition <- matrix(mTransition,nrow=5,ncol=5);
mTransition <- t(mTransition)
## 2. Steady State
capitalSteadyState = (aalpha*bbeta)^(1/(1-aalpha));
outputSteadyState = capitalSteadyState^aalpha;
consumptionSteadyState = outputSteadyState-capitalSteadyState;
cat(" Output = ", outputSteadyState,", Capital = ",capitalSteadyState, ", Consumption = ", consumptionSteadyState,"\n")
cat(" \n")
# We generate the grid of capital
vGridCapital <- seq(0.5*capitalSteadyState, 1.5*capitalSteadyState, by = 0.00001);
nGridCapital <- length(vGridCapital);
nGridProductivity <- length(vProductivity);
# 3. Required matrices and vectors
mOutput <- matrix(0,nGridCapital,nGridProductivity);
mValueFunction <- matrix(0,nGridCapital,nGridProductivity);
mValueFunctionNew <- matrix(0,nGridCapital,nGridProductivity);
mPolicyFunction <- matrix(0,nGridCapital,nGridProductivity);
expectedValueFunction <- matrix(0,nGridCapital,nGridProductivity);
## 4. We pre-build output for each point in the grid
mOutput = as.matrix(vGridCapital^aalpha)%*%t(as.matrix(vProductivity));
## 5. Main iteration
maxDifference <- 10;
tolerance <- 0.0000001;
iteration <- 0;
while (maxDifference>tolerance){
expectedValueFunction = mValueFunction %*% t(mTransition);
for (nProductivity in 1:nGridProductivity){
# We start from previous choice (monotonicity of policy function)
gridCapitalNextPeriod <- 1;
for (nCapital in 1:nGridCapital){
valueHighSoFar <- -100000;
capitalChoice <- vGridCapital[1];
for (nCapitalNextPeriod in gridCapitalNextPeriod:nGridCapital){
consumption <- mOutput[nCapital,nProductivity]-vGridCapital[nCapitalNextPeriod];
valueProvisional <- (1-bbeta)*log(consumption)+bbeta*expectedValueFunction[nCapitalNextPeriod,nProductivity];
if (valueProvisional>valueHighSoFar){
valueHighSoFar <- valueProvisional;
capitalChoice <- vGridCapital[nCapitalNextPeriod];
gridCapitalNextPeriod <- nCapitalNextPeriod;}
else{
break; # We break when we have achieved the max
}
}
mValueFunctionNew[nCapital,nProductivity] <- valueHighSoFar;
mPolicyFunction[nCapital,nProductivity] <- capitalChoice;
}
}
maxDifference <- max(abs(mValueFunctionNew-mValueFunction));
mValueFunction <- mValueFunctionNew;
iteration = iteration+1;
if ((iteration %% 10)==0 | iteration ==1){
cat(" Iteration = ", iteration," Sup Diff = ", maxDifference,"\n");
}
}
cat(" Iteration = ", iteration," Sup Diff = ", maxDifference,"\n");
cat(" \n")
cat(" My chek = ", mPolicyFunction[1000,3],"\n");
cat(" \n")
cat(" Time = ", proc.time() - ptm,"\n");