# "Array indices must be positive integers or logical values." while using golden section search with an interpolation?

3 visualizaciones (últimos 30 días)
Kevin Zhou el 1 de Nov. de 2020
Comentada: Kevin Zhou el 1 de Nov. de 2020
The function starts going wrong as I hit line 50, the one that starts with w(i,j). What am i doing wrong? A seperate function is implemented in the code called gss which represents the golden search function, code posted below:
%This script performs value function iteration for the baseline model.
clear all
global delta eta alpha beta
delta=1; %depreciation rate
eta=1; %utility parameter
alpha=.27; %production function parameter
beta=.994; %patience parameter
eps=.01; %partial tolerance parameter
n=250; %grid size
Kl=.75*kstar; %lower bound for capital grid
Ku=1.25*kstar; %upper bound for capital grid
kstep=(Ku-Kl)/(n-1); %step size for capital grid
kgrid=Kl:kstep:Ku; %grid for capital stock
v0=zeros(n,1); %initialize the value function
dist=1; %initialize the distance for the value function
it=1; %initialize the counter
tic %initialize the clock
eps2=.000001;
while (dist>eps*(1-beta))
for i=1:length(kgrid)
%Step 1
for j=1:length(kgrid)
ct=kgrid(i)^alpha+(1-delta)*kgrid(i)-kgrid(j); %consumption
%if i=1, kgrid(i)=0.1237. If I choose a kgrid of j=1, the next
%period will also be 0.1237. As we mvoe to the next period,
%kgrid(j)=0.1246.
ct=ct*(ct>0)+eps2; %correct for potential negative consumption
%ij=[kgrid(j-1) krid(j+1)]- needed component;
%yi_given=interpl(x,y,xi_given,'linear'); x&y are given strings
%of values or you can represent it as below
%yp=interpl(x,y,xp,'linear'); xp=[0:0.25:10]; The interpolation
%is a way to make the points finer to match the now finer gird
%neeeded for golden section searhc as the gss method demand
%the input of an interval with only upper and lower bound with
%no step size.
%you set i=1 and go through all these js. For every single i,
%you will try to find the j that maxmizes that value function.
%Essentially, you keep the i fixed, i=1 and go through all the
%columns represneted by j. For instance, it goes from (1,1) to
% (1,2 & n) for i=1. You are trying to find the max within that
%row through step 2.1.
w(i,j)=gss(utilchp4(ct,eta)+beta*interpl(1:j,v0(j,1),kgrid(j-1):krid(j+1),'linear'),[kgrid(j-1) krid(j+1)]);
%value function are
%monotonically increasing. This function stores the value
%created by capital as indicated by i and j.utilchp4(ct,eta)
%give you a instaneous return. Use golden section search to
%lienar interpolate that entire right hand side. Professor
%highlighted v0(j,1) when speaking about linear interpolation.
end
%Step 2
%find the max of w for every i. When you find that, that points to
%jstar. The location of the column of the biggest W essentially
%represent the location of the column of the biggest j.
jind=find(max(w(i,:))==w(i,:)); %find the j that max rhs of Bellman.
%when i=2, jstar will be 1, which indicate the first value in the
%whole string of values as laid out by kgrid.
vnew(i,1)=w(i,jind); %update the value function for each i. This
%tells you where in W that max is.
%Step 3
h(i,1)=jind; % form your policy function. For every i, store jstar.
%For each iteration, store jstar. Check line 85.
end
dist=max(abs(v0-vnew)); %calculate distance between value functions
diststore(it)=dist; %store the distance between value function iterations
%pause(.2)
%plot(1:1:it,diststore)
it=it+1; %update counter
v0=vnew;%update the value function and repeat. The method discussed in
%in class notes as we keep on using the last kprime to get converged
%toward the kprime
end
kprime=kgrid(h); %given policy pointer h, find the capital stock value
disp('time to completion')
toc
The golden section search code is provided below (validated with no problem):
function [x]=gss(fun,X)
eps= 0.000001;
p=(sqrt(5)-1)/2; %the golden ratio
q=1-p;
% Step 1
a=X(1); %original lower bound
d=X(length(X)); %original upper bound
b=p*a+q*d; %lower bound combined with golden ratio calculation
c=q*a+p*d;%upper bound combined with golden ratio calculation
fb=fun(b);
fc=fun(c);
% step2+3
while abs(d-a)>eps*max(1,(abs(a)+abs(c))) %make the interval smaller and smallers
if fb<fc % choose [b,d] as next interval
a=b; %cut the range from the right hand side
b=c;
fb=fc;
c=p*b+q*d;
fc=fun(c);
else % choose [a,c] as next interval
d=c; %cut the range from the left hand side
c=b;
fc=fb;
b=p*c + q*a;
fb=fun(b);
end
end
%as the function keeps on iterating....
if fb>fc
x=b;
else
x=c;
end
%disp('To prove that it actually equals to 0')
%the_function_result=fun(x)
end
##### 0 comentariosMostrar -2 comentarios más antiguosOcultar -2 comentarios más antiguos

Iniciar sesión para comentar.

### Respuestas (1)

Walter Roberson el 1 de Nov. de 2020
w(i,j)=gss(utilchp4(ct,eta)+beta*interpl(1:j,v0(j,1),kgrid(j-1):krid(j+1),'linear'),[kgrid(j-1) krid(j+1)]);
What are you expecting the kgrid portion of that code to do when j is 1 or j is length(kgrid)?
##### 1 comentarioMostrar -1 comentarios más antiguosOcultar -1 comentarios más antiguos
Kevin Zhou el 1 de Nov. de 2020
That is a good question. Idealistically, I just want them to start and stop at j is 1 and lenght(kgrid). Would you know how to modify it?

Iniciar sesión para comentar.

### Categorías

Más información sobre Matrix Indexing en Help Center y File Exchange.

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by