Life-cycle models with iid shock

When using iid shocks, you have to pass the grid e_grid and probability vector pi_e to vfoptions and simoptions.

I have passed pi_e as a row vector to the distribution by mistake, i.e.

vfoptions.pi_e   = pi_e(1,:)';
simoptions.pi_e   = pi_e(1,:); % Here I forgot the transpose!!!

and the code returns with an error

Error in StationaryDist_FHorz_Case1 (line 161)
simoptions.pi_e_J=simoptions.pi_e.*ones(1,N_j,‘gpuArray’);

This can be easily fixed by typing

simoptions.pi_e   = pi_e(1,:)'

Interestingly, if I pass a row vector pi_e to the VFI, then I get an informative error:

vfoptions.pi_e is not the correct shape (should be of size N_e-by-1 or N_e-by-N_j)

I would propose to add this error also to the stationary distribution, or to add to both VFI and distribution the following three lines:

if isrow(simoptions.pi_e)
   simoptions.pi_e=simoptions.pi_e'
end

Unrelated to the error: Is it possible to pass two iid shocks? My idea is

vfoptions.e_grid = [e1_grid;e2_grid];
vfoptions.pi_e = kron(e2_prob,e1_prob); %note kron in reverse order
1 Like

The reason it doesn’t check for this is I figured everyone would just use:

simoptions.e_grid=vfoptions.e_grid; 
simoptions.pi_e=vfoptions.pi_e;

and hence as long as I checked and warned about the shape in value fn, no-one would get it wrong in agent dist so I didn’t need to add extra overhead checking it again (you cannot run agent dist till after you have policy from running value fn, and to use simoptions.e_grid you must also be using vfoptions.e_grid).

Clearly though this is not how your code was written. Is there a reason for you not just using simoptions.pi_e=vfoptions.pi_e; If there is a reason I will add the check to agent dist codes?, but otherwise I will probably just leave it as is (otherwise it adds more overhead).

Two shocks:

vfoptions.e_grid=[e1_grid; e2_grid];
vfoptions.pi_e=kron(pi_e2, pi_e1); % note reverse order
simoptions.e_grid=vfoptions.e_grid;
simoptions.pi_e=vfoptions.pi_e;

PS. The other nice thing is that you can just mindlessly copy-paste
simoptions.e_grid=vfoptions.e_grid;
simoptions.pi_e=vfoptions.pi_e;
into every code where you use e shocks :smiley:

1 Like

Actually, as a longer term goal, I really want to make it so that lots of these ‘checks’ can be skipped in things like general eqm commands where you are calling these commands lots of times. If this were the case then I would not be worried about the overhead from things like checking e_grid in simoptions after it has already been done in vfoptions.

1 Like

No reason in particular, other than a useless repetition :rofl:
It makes sense not to clutter the code with too many checks

The line

vfoptions.pi_e=kron(e1_grid, e2_grid);

Should read as

vfoptions.pi_e=kron(e1_prob, e2_prob);

where e1_prob and e2_prob are column vectors with the probability weights of e1 and e2, respectively. Thanks a lot for point out the kron ordering, I had completely overlooked it!

1 Like

fixed my typo, thanks.

1 Like

I’m not sure I am convinced about the ordering of the kron :slight_smile:

Say you have x=(x1,x2) and y=(y1,y2) as iid shocks, with probabilities p_x=(p_x1,p_x2) for x and p_y=(p_y1,p_y2) for y. In general, there are nx points for x and ny points for y (here for simplicity, nx=ny=2).

I choose to vary x first, hence the combined variable e (with size [nx*ny,2]) will be as follows:

e1 = x1,y1
e2 = x2,y1
e3 = x1,y2
e4 = x2,y2

The corresponding probability vector p_e (with size [nx*ny,1]) will be as follows:

p_e1 = p_x1 * p_y1
p_e2 = p_x2 * p_y1
p_e3 = p_x1 * p_y2
p_e4 = p_x2 * p_y2

Using the kron product, I can write p_e as a function of p_x (size [nx,1]) and p_y (size [ny,1]) as follows:

p_e = p_y \kron p_x = kron(p_y,p_x)

because

p_y \kron p_x = [ p_y1*[p_x1,
                                         p_x2]
                            p_y2*[p_x1,
                                        p_x2] ]

which gives a 4*1 column vector, with the p_x elements changing first. The ordering is e=(x,y) but the kron is in reverse order, as in the standard case. Whether it is row or column should not matter.

The definition of kron product is here: kron

Please let me know what you think, especially in the likely case that I missed something!

2 Likes

Sorry, you were right, kron() should have been reverse order. Fixed my post.

Note sure what happened there, I had created aa=rand(3,1) with bb=[2;1] and then looked at kron(aa,bb) and kron(bb,aa). Somehow my brain shut off and I chose the wrong one! Anyway, fortunately you were awake :slight_smile: Sorry for wasting a bit of your time there by giving you the wrong directions :pray:

Actually nowadays I often avoid kron() because I don’t find it intuitive, instead I tend to use repmat() and repelem() which I find much more intuitive and which can be used in high-dimensional setups. For example I could create pi_e as pi_e=repmat(pi_e1,n2,1).*repelem(pi_e2,n1,1);. More generally you just repmat() for any tailing parts in that dimension and repelem() for any preceding parts in that dimension, you can use both nested if you have both preceding and tailing parts. That said, maybe I just find these easier because I am doing lots of high dimensional stuff inside the toolkit, and for a standard user it is probably just easier to remember to use kron() with reverse ordering to combine transition matrices.

I copy paste something from your post here, because it gives a nice clean explanation of what we want to happen


I choose to vary x first, hence the combined variable e (with size [nx*ny,2]) will be as follows:

e1 = x1,y1
e2 = x2,y1
e3 = x1,y2
e4 = x2,y2

The corresponding probability vector p_e (with size [nx*ny,1]) will be as follows:

p_e1 = p_x1 * p_y1
p_e2 = p_x2 * p_y1
p_e3 = p_x1 * p_y2
p_e4 = p_x2 * p_y2
2 Likes

No problem, thanks for the clarification! I didn’t know the repelem function, thanks for suggesting it. Indeed, kron works only with vectors and matrices

For multidimensional kron product, there is a function written by Mathias Kredler: