Making interpolations in a 4D matrix

I have Matlab code that is very inefficient and I need to run it several times. The code is basically a big parfor loop which I guess there is almost impossible to get around.

The code starts by loading several parameters and 4-D matrices, and then needs to make a few interpolations. All need to be done 5000 times (thus the parfor loop).

I simplified the code the most I could without taking out the key ingredients.

load file

nsim = 5000
T = 12;
N = 1000;

cumQx = cumsum(Qx);
cumQz = cumsum(Qz);
cumQs = cumsum(Qs);

for k=1:nsim
st(k).ksim    = kstar*ones(N, T);
st(k).Vsim  = zeros(N,T);
st(k).Psim = zeros(N,T);
end

parfor k = 1:nsim

sysrand  = rand(T, 1);
idiorand = rand(N, T);
sigmarand = rand(T,1);

xid = zeros(T, 1);
zid = zeros(N, T);
sid = zeros(T,1);
xid(1) = 8;
zid(:, 1) = 5;
sid(1) = 1;
% Initializing the simulation

simx    = zeros(T,1);
zsim    = ones(N,T)*zbar;
simsx    = zeros(T,1);

% Construct 3-D grid using 'ndgrid'
[ks,zs] = ndgrid(kgrid,z);

for j = 2:T
sid(j) = find(cumQs(:, sid(j-1)) >= sigmarand(j), 1);
simsx(j-1) = sigmax(sid(j));

xid(j) = find(cumQx(:, xid(j-1)) >= sysrand(j), 1);
simx(j-1) = x(xid(j));

for n = 1:N
zid(n, j)   = find(cumQz(:, zid(n, j-1)) >= idiorand(n, j), 1);
zsim(n,j-1) = z(zid(n, j));
end
st(k).ksim(:,j)      = interpn(ks, zs , squeeze(kprime(:,xid(j),:,sid(j))),   st(k).ksim(:,j-1),zsim(:,j-1),'linear');       % K
st(k).Vsim(:,j)      = interpn(ks, zs , squeeze(V(:,xid(j),:,sid(j))),        st(k).ksim(:,j-1),zsim(:,j-1),'linear');       % V
st(k).Psim(:,j)      = interpn(ks, zs , squeeze(P(:,xid(j),:,sid(j))),        st(k).ksim(:,j-1),zsim(:,j-1),'linear');       % P

end;

end


Here is a link for the matrices needed to run the code. Also a dropbox link here.

Is there any better way of doing this that would significantly reduce the computation time? My guess is no.