2-d interpolation in matlab

I have a matrix v0 nd*nk*na, which depends on values of d, k and a. I want to interpolate the matrix v for values of d and k which are stored in dvec and kvec and range from dmin to dmax and kmin and kmax respectively. I tried to use interp2d as follows:

g = interp2(kvec,dvec,v0,k,d,'linear');

I get this error:

Error using  .' 
TRANSPOSE does not support N-D arrays. Use PAGETRANSPOSE/PAGECTRANSPOSE to transpose pages or PERMUTE to reorder
dimensions of N-D arrays.

Is there any way to fix this?

[Matt J]

You haven't shown the complete error message (why not?), but it doesn't appear to be coming from interp2.

[Prerna Mishra]

This is the next part:

Error in interp2 (line 120)
V = V.';

[John D'Errico]

It appears you are using interp2 in a way that it is not designed to be used. Interp2 allows you to index into a 2 dimensional matrix essentially at intermediate points between elements.

However, v0 is a THREE dimensional matrix. So you would use THREE variables to index into it, yet you want to provide only TWO variables. Do you see the problem?

[Prerna Mishra]

I see. But is there any way to do the problem that I have? Keeping the a dimension constant, can I interpolate the value of v0 for d and k?

[John D'Errico]

@dpb has given you the correct answer, where I was being too lazy to go. ;-)

Effectively, you need to treat this as a 3-dimensional interpolation. Or you can extract each plane of v0, and then use a 2-d interpolation along that plane, thus for each possible a.

Answer 1

[dpb]

As per usual, would be far easier to visualize what you're after if would provide a (smallish) sample dataset.

But, if by " Keeping the a dimension constant, ..." you mean interpolating over the other two dimension at a given plane of the 3D array, then sure -- just pass the desired plane into the interp2d call...

ixa=...       % a specific plane of 3D array v0 in range 1:size(v0,3)
g = interp2(kvec,dvec,v0(:,:,ixa),k,d,'linear');

Alternatively, coding may be simpler to just use interp3 even if the 3rd dimension interpolated value is constant.

g=interp3(kvec,dvec,avec,v0,k,d,a,'linear');

Nothing says a can't be a single value.

[Prerna Mishra]

I tried this:

g=interp3(dvec,kvec,avec,v0,d,k,a0,'linear'); %v0 is nd*nk*na, 
% d is in the first dimension, k in the second and a in the third. 

a0 is the constant.

I now get the following error:

Error using griddedInterpolant
Sample points vector corresponding to grid dimension 1 must contain 11 elements.

Error in interp3 (line 132)
        F = griddedInterpolant({X, Y, Z}, V, method,extrap);

[Torsten]

Please run the following code and tell us what MATLAB gives as output concerning the sizes:

size(dvec)
size(kvec)
size(avec)
size(v0)
size(d)
size(k)
size(a0)
g=interp3(dvec,kvec,avec,v0,d,k,a0,'linear');

[Prerna Mishra]

I was able to fix this, sorry just a size mismatch.

[dpb]

Hmmm....I don't recall I had ever used interp3 without building a meshgrid [X,Y,Z] array, but in keeping with the cartesian coordinate orientation, interp3 is expecting X to be the horizontal and Y to be vertical; hence use

g=interp3(kvec,dvec,avec,v0,d,k,a0,'linear');

instead.

Illustrating with a sample case

M=rand(6,7,3);
X=1:6;Y=1:7;Z=1:3;
interp3(X,Y,Z,M,1,1,1)

Error using griddedInterpolant
Sample points vector corresponding to grid dimension 1 must contain 7 elements.

Error in interp3 (line 132)
F = griddedInterpolant({X, Y, Z}, V, method,extrap);

So, need to use

interp3(Y,X,Z,M,1,1,1)

instead

[Prerna Mishra]

I get a another error. For simplification, I have written ta simplfied version of the entire code here. Basically it is a stochastic dynamic programming problem, trying to optimize over levels of capital and debt. The main function uses fminsearchbnd to find the value of k and d that minimizes the value function. The valfun_stoch uses the interpolation since fminsearchbnd will not necessarily find the values of capital and debt that I defined in the grid kvec and dvec. Following is the code I am running,

main:

clear all;
close all;
tic
global v0 beta delta alpha kmat k0 s prob a0 s j val g v1 i dmat amat kd

%set up the fised paramters
alpha = 0.33; %capital share
beta = 0.95;%discount rate
delta = 0.1; %rate of depreciation of capital
s = 2;


tol = 0.01;%tolerance
maxits = 1000;%maximum iteration
dif = tol + 100;
its = 0;

kstar = (alpha/(1/beta-(1-delta)))^(1/(1-alpha));
cstar = kstar^(alpha)-delta*kstar;
istar = delta * kstar;
ystar = kstar ^ alpha;

kmin = 0.25*kstar;
kmax = 1.75*kstar;
kgrid = 99;
grid = (kmax - kmin)/kgrid;
kmat = kmin:grid:kmax;
kmat = kmat';
dmin = kmin/2;
dmax = kmax/2;
dgrid = 99;
d_grid = (dmax - dmin)/dgrid;
dmat = dmin:d_grid:dmax;
dmat = dmat';
[N,n] = size(kmat);

%polfun = zeros(kgrid+1,3);
v0 = zeros(N,N,3);

amat = [0.9 1 1.1]';
prob = (1/3) * ones (3,3);


while dif>tol && its < maxits
    
    
    for j = 1:3
        for i = 1:N
            for l = 1:N 
                k0 = kmat(i,1);
                a0 = amat(j,1);
                
                kd = fminsearchbnd(@valfun_stoch,[kmin dmin],[kmin dmin],[kmax dmax]);
                %fun = @(k) valfun_stoch;
                
                %k1 = fminsearch(@valfun_stoch,kset);
                %k1 = fminbnd(@valfun_stoch,kmin,kmax);
                v1(l,i,j) = -valfun_stoch(kd);
                k11(l,i,j) = kd(1);
            end 
        end
    end
    %g = abs(v1−v0);
    dif = norm(v1-v0);
    v0 = v1;
    its = its+1;
end
toc

and the valfun_stoch is as follows:

function val=valfun_stoch(kd)

global v0 beta delta alpha kmat k0 prob a0 s j val g v1 i dmat amat

g = interp3(kmat,dmat,amat,v0,kd(1),kd(2),a0,'linear'); % smooths out previous value function

c = a0*k0^alpha - kd(1) + (1-delta)*k0 - kd(2); % consumption
if c <= 0
    val = -8888888888888888 - 800*abs(c); % keeps it from going negative
else
    val = (1/(1 - s))*(c^(1 - s) - 1) + beta*(g*prob(j,:)');
    
end
val = -val; %

The error is:

Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 3-by-1.

Error in fminsearch (line 201)
fv(:,1) = funfcn(x,varargin{:});

Although at this moment it is quite possible that fminsearchbnd is throwing an error and not interpolation

[dpb]

Without the full error in context we really can't tell because of the excessive use of global where the size mismatch actually is in your user code; the error message you posted is internal.

What it's telling you though, is that somewhere in your code is something like

...
results=[1:3].';        % internal code generates a three-vector somehow
y(1)=results;  % user code tries to stuff a vector into a single array element -- BOOM!

If I had to guess, I'd say it's probably

v1(l,i,j) = -valfun_stoch(kd);

where valfun is returning a three vector.

Use the debugger to find your logic/coding error in not maintaining dimensions...use either

dbstop on error

or just set a breakpoint inside valfun and step through it to be sure it does what is expected.

[Torsten]

@Prerna Mishra

Well, I told you several times now that the assignment

val = (1/(1 - s))*(c^(1 - s) - 1) + beta*(g*prob(j,:)');

makes "val" a (3x1) vector. But "val" must be a scalar (1x1) for the optimizer to work.

[Prerna Mishra]

Right, I understood that. The problem is that when I interpolate using the code above, I get g which is a scalar. I should have ideally gotten a 7 * 1 vector.

For example in an easier problem say:

v0 = rand(100,3);
kv = 1:100;
g = interp1(kv,v0,1.5,'linear')

The output of interpolation is 3*1.

I was expecting similar for this 3 - dimensional problem. That the output of the interpolation should be na*1 vector and not a scalar.

[Torsten]

Let the result of the interpolation be whatever it is, but in the end, "val" must be a scalar.

If g is a vector of the same size as prob(j,:), g*prob(j,:).' gives a scalar - so this would work.

Since we all don't know what you are doing in your code, we can only point out the mistakes, but cannot give advice on how to change your code adequately.