Calculate Correlations (and Covariances) from values, weights, and transition probabilities

1

Following is Matlab code demonstrating how to compute autocorrelation coefficients from matrices of values and weights, as well as the transition probabilities [so for a discrete markov process].


% Calculate autocorrelations from values, weights, and transition probabilities

% Idea: there are N points in state-space
% We have weights W and values V (on the state-space)
% We have transition probabilities P
% [In finite horizon, or just generally outside the a stationary eqm, we
% also have next period weights W2 and values V2. In a stationary eqm for
% infinite horizon model these will just be copies of W and V.]
%
% First section of code just sets all this up
% Second section of code shows how to calculate the correlation of the two
% periods. Done by first calculating the covariance and the standard
% deviations.

%% Setup
InfHorz=1;
HighCorr=0; % change P to make correlation high (only works for infinite horizon setup)

% Number of points in state space
N=2;

% We have function of the state space, and this period it takes values
V=10*rand(N,1);

% Transition matrix for the states is
P=rand(N,N);
P=P./(sum(P,2)); % normalize rows to 1
if HighCorr==1 % Force high correlation
    P=P+eye(N,N); % move mass towards diagonal
    P=P./(sum(P,2)); % normalize rows to 1
elseif HighCorr==2 % perfect correlation, so will get corr=1
    P=eye(N,N);
end


% Set up weights for this period and next period, and values for next period
if InfHorz==1 % Infinite horizon
    % Use the weights that correspond to the stationary dist of P
    dist=1;
    W=rand(N,1);
    W=W/sum(W); % total mass of 1
    while dist>0
        Wlag=W;
        W=P'*W;
        dist=max(abs(W-Wlag));
    end

    V2=V;
    W2=W;

    fprintf('double-check, following should be zero (or at least 10^-9 ish): %1.10f \n', max(abs(W-(P'*W))) )
else
    V2=10*rand(N,1);

    % Say the distribution over states this period is
    W=rand(N,1);
    W=W/sum(W); % total mass of 1
    % Then value over states next period is
    W2=P'*W;

end



%% Autocorrelation
% We want to calculate the autocorrelation
% Correlation(x,y)=Cov(x,u)/(stddev(x)*stddev(y))
% So first calculate the covariance and the two standard deviations

meanV=sum(W.*V);
meanV2=sum(W2.*V2);
stddevV=sqrt(sum(W.*(V-meanV).^2));
stddevV2=sqrt(sum(W2.*(V2-meanV2).^2));

% Calculate covariance between this period and next period values
Covar=(W.*V)'*P*V2 - meanV*meanV2;
fprintf('Covariance is: %2.4f  \n', Covar)

% Calculate the correlation
Corr=Covar/(stddevV*stddevV2);
fprintf('Correlation is: %2.4f  \n', Corr)


%% Double check by simulating the model and calculating the autocorrelation from this
S=10000; % number of simulations

% first, get pairs of indices for this period and next period
Wcum=cumsum(W);
Pcum=cumsum(P,2);
indexThisPeriod=zeros(S,1);
indexNextPeriod=zeros(S,1);
for ss=1:S
    draw=rand(1,1);
    [~,ii]=max(Wcum>draw,[],1);
    indexThisPeriod(ss)=ii;
    draw2=rand(1,1);
    [~,ii2]=max(Pcum(ii,:)>draw2,[],2);
    indexNextPeriod(ss)=ii2;
end
Vsim=V(indexThisPeriod);
V2sim=V2(indexNextPeriod);

CovarSim=cov(Vsim,V2sim);
CovarSim=CovarSim(1,2); % the covariance is the off-diagonal [Matlab computes the whole Covariance Matrix, the diagonals are the variances]

fprintf('Covariance from two different approaches should be roughly same: %2.4f, %2.4f  \n', Covar, CovarSim)
% If you set S really large, they should become essentially identical.
% Only difference should be the sampling error


% Double-check: W and histcounts(indexThisPeriod) should be roughly the same (once you normalize the later; looks fine)
% H=histcounts(indexThisPeriod)'; H/sum(H)
% Double-check: if you set HighCorr==2 then get indexThisPeriod=indexNextPeriod (passed this)
1 Like
2

If you want, e.g., five-period autocorrelations you just use P^5 in place of P (the five-period transition matrix in place of the one-period transition matrix), or the finite-horizon equivalent of doing this.

3

This is really useful. Going back to @aledinola example, I can set

W=W2=StationaryDist(:slight_smile:
V=V2=pol_n(:slight_smile:

How do I define P? Is it simply pi_z?

1 Like
4

Hi @javier_fernan,

I think that P is Robert’s code is a matrix that gives the probability of going from x to x' where x is a generic point on the state space. In the example, x=(a,z). So you can think of P(x,x')=Prob(a',z' | a,z). To generate this matrix you have to combine the policy function for a' with the exogenous transition \pi(z,z').
Here is my code to do this:

Ptemp = cell(N_z,1);
for z_c=1:N_z
   Ptemp{z_c} = sparse((1:N_a),Policy(:,z_c),1,N_a,N_a);
   Ptemp{z_c} = kron(pi_z(z_c,:), Ptemp{z_c} );
end
P = cell2mat(Ptemp);

You can verify that (1) the size of P is [N_a*N_z,N_a*N_z] and that (2) the rows of P sum to one.

1 Like
5

As @aledinola says, P will map (a,z) into (a’,z’), and is constructed by combining Policy with pi_z.

[Note: this is the same P you would use to iterate on the agent distribution. Of course, nowadays you would use Tan improvement instead to do the iterations as two sub-steps without actually constructing the whole P. But just mention in case you ever saw that before.]

1 Like