function [invm,invmt,vt,vr,copt,colr,sopt,solr,avgk,aveinc]=...
    vfi_retirement_TL_transitory_pens(dar,dpol,pens,ai,ar,sig,...
                bet,intpcf,ylife,Qylife,r,maxa,...
                invherhip,ye,yedtb,initialwealth,...
                blimitage)
            
 global da dy dbet TL T TR gpop a sur dhip dye dagg azero
           
%%%%%%%%%%%%%%%% value and policy function initializations %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% retired agents
vr=zeros(da,dy*dye,dbet,TR+1);

% initialisation of the retired agent's policy function 
polr=zeros(da,dy*dye,dbet,TR);
colr=polr;

%%%%%%%%%%%%%%%% working people
vt=zeros(da,dye,dy,dhip,dbet,TL);
copt =vt; % optimal consumption policy function
sopt=vt; % optimal savings policy function.

% This is now disposable income, no taxation

    ttau=0;
    tlamb=1;
    beq=0; % bequests are thrown to the ocean

   %Now start solving backwards for the value functions
   
%% 6.2. Retired agents

for j=TR:-1:1
   dc=100;
      
      for jjj=1:dy*dye   % jjj is the index for y_t
          
          
      vri=[vr(:,jjj,:,j+1); intpcf*vr(da-dpol+1:da,jjj,:,j+1)];

         
         % first compute the utility matrix for each value of consumption
         % (rows)and today's assets (columns)
         c=linspace(0,1,dc+1)'; % proportion of this period's income consumed
         %c=c(2:dc+1)*((1+r*(1-taua))*a'+(1-taul)*(...
         %    ylife(dy*dhip*(TL-1)+dy*(j2-1)+jjj))+beq); % actual
         % consumption values
         
         % Notice that now you can consume more than 100%
         
         %c=c(2:dc+1)*(a'+tlamb*(pens(jjj)+r*a').^(1-ttau)+beq);
         c=c(2:dc+1)*(a'+blimitage(TL+j)+tlamb*(pens(jjj)+r*a').^(1-ttau)+beq);
         c(c<0)=NaN;

         u=c.^(1-sig)/(1-sig);

         % compute tomorrow's assets when not inheriting (=savings) for each
         % possible value of consumption and assets today
         s=-c+ones(dc,1)*(a'+blimitage(TL+j)+tlamb*(pens(jjj)+r*a').^(1-ttau)+beq)-...
             blimitage(TL+j);
         % reshape them in order to interpolate the value function later on
         st=reshape(s,1,da*dc);

         % compute tomorrow's assets when inheriting, for each possible level
         % of bequests received tomorrow (on the same grid as assets)
         % and each possible level of consumption
         %at=a*ones(1,da*dc)/numch+ones(da,1)*st;
         % reshape them in order to interpolate the value function
         %at=reshape(at,1,da^2*dc);

         for jb=1:dbet % j is the index for beta

            %%% compute the value function when inheriting
            % compute the value function for all the possible levels
            % of assets tomorrow (at)
            %vat=interp1q(ai,vri(:,j),at')';
            % now reshape the interpolated value function to compute its
            % expected value with respect to the bequests later on
            %vat=reshape(vat,da,da*dc);

            %%%% compute the value function when not inheriting
            vst=interp1q(ai,vri(:,jb),st')';
            % reshape again
            vst=reshape(vst,dc,da);
      

            [temp1, temp2]=...
            max(u+sur(TL+j)*bet(jb)*vst);
            vr(:,jjj,jb,j)=temp1';
            polr(:,jjj,jb,j)=temp2';
            

            
   
            colr(:,jjj,jb,j)=(temp2./dc)'.*(a+blimitage(TL+j)+...
               tlamb*(pens(jjj)+r*a).^(1-ttau)+beq);
                              %temp2 refers to choice of
                              %today's CONSUMPTION (it is what we are
                              %substituting out here). We divide by dc in
                              %order to find out the proportion of today's
                              %assets that is being consumed and then we
                              %multiply by today's assets, finding finally
                              %the level of consumption.
                              
                              %Note c was 0:dc+1 but then we just picked
                              %c(2:dc+1). This leaves a c where proportions
                              %are (for dc=100) 0.01, 0.02... 1, which
                              %means temp2/dc picks the right proportion as
                              %well.
                  
           
            solr(:,jjj,jb,j)=(ones(da,1)-(temp2./dc)').*...
                                  (a+blimitage(TL+j)+...
               tlamb*(pens(jjj)+r*a).^(1-ttau)+beq)-blimitage(TL+j);
           
        if (colr(:,jjj,jb,j)<0)
            colr(:,jjj,jb,j)=NaN;
            vr(:,jjj,jb,j)=-1e5;
        end
        
       vr(isnan(vr(:,jjj,jb,j)),jjj,jb,j)=-1e5;

                
         end % for jb=1:dbet, j is for bet
       end % for jjj=1:dy, jjj is for h_t
 
end
   

   %% 6.3. Working agents

   % we need a grid for consumption, since we are subtituting a_{t+1}
   % in this case; the range for consumption will be different for different
   % asset level, productivity level and age
   % Grid of consumption is the proportion of assets consumed in a period
   % define consumption grid size
   dc=100;

   %%%%%%%%%%%% age TL is special: next period the guy is retired and the
      % corresponding value function has only 3 dimensions
      
      % let's take care of a "big" vr for the agent who will be retired
      % next period (in order not to fall out of the interpolation range)
     for j2=1:dhip %j2 is the index for the hip group
      for jjj=1:dy   % jjj is the index for y_t
      for je=1:dye %je is the index for eps_t
          
          % Note here the order of pens (and of dy*dye) matters
          % Let's consider the transitory shock varies faster

          
      vri=[vr(:,(jjj-1)*dye+je,:,1); intpcf*vr(da-dpol+1:da,(jjj-1)*dye+je,:,1)];

         
         
         % first compute the utility matrix for each value of consumption
         % (rows)and today's assets (columns)
         c=linspace(0,1,dc+1)'; % proportion of this period's income consumed
         %c=c(2:dc+1)*((1+r*(1-taua))*a'+(1-taul)*(...
         %    ylife(dy*dhip*(TL-1)+dy*(j2-1)+jjj))+beq); % actual
         % consumption values
         c=c(2:dc+1)*(a'+blimitage(TL)+tlamb*(ylife(dy*dhip*(TL-1)+...
             dy*(j2-1)+jjj)*ye(je,TL)+r*a').^(1-ttau)+beq);
         
         c(c<0)=NaN;

         u=c.^(1-sig)/(1-sig);

         % compute tomorrow's assets when not inheriting (=savings) for each
         % possible value of consumption and assets today
         s=-c+ones(dc,1)*(a'+blimitage(TL)+tlamb*(...
             ylife(dy*dhip*(TL-1)+dy*(j2-1)+jjj)*ye(je,TL)+r*a').^(1-ttau)+beq)-...
             blimitage(TL);
         % reshape them in order to interpolate the value function later on
         st=reshape(s,1,da*dc);

         % compute tomorrow's assets when inheriting, for each possible level
         % of bequests received tomorrow (on the same grid as assets)
         % and each possible level of consumption
         %at=a*ones(1,da*dc)/numch+ones(da,1)*st;
         % reshape them in order to interpolate the value function
         %at=reshape(at,1,da^2*dc);

         for j=1:dbet % j is the index for beta

            %%% compute the value function when inheriting
            % compute the value function for all the possible levels
            % of assets tomorrow (at)
            %vat=interp1q(ai,vri(:,j),at')';
            % now reshape the interpolated value function to compute its
            % expected value with respect to the bequests later on
            %vat=reshape(vat,da,da*dc);

            %%%% compute the value function when not inheriting
            vst=interp1q(ai,vri(:,j),st')';
            % reshape again
            vst=reshape(vst,dc,da);
      

            [temp1, temp2]=...
            max(u+sur(TL)*bet(j)*vst);
            vt(:,je,jjj,j2,j,TL)=temp1';
            %polt(:,je,jjj,j2,j,TL)=temp2';
            
            vt(isnan(vt(:,je,jjj,j2,j,TL)),je,jjj,j2,j,TL)=-1e5;
            
   
            copt(:,je,jjj,j2,j,TL)=(temp2./dc)'.*(a+blimitage(TL)+...
               tlamb*(ylife(dy*dhip*(TL-1)+dy*(j2-1)+jjj)*ye(je,TL)+r*a).^(1-ttau)+beq);
                              %temp2 refers to choice of
                              %today's CONSUMPTION (it is what we are
                              %substituting out here). We divide by dc in
                              %order to find out the proportion of today's
                              %assets that is being consumed and then we
                              %multiply by today's assets, finding finally
                              %the level of consumption.
                              
                              %Note c was 0:dc+1 but then we just picked
                              %c(2:dc+1). This leaves a c where proportions
                              %are (for dc=100) 0.01, 0.02... 1, which
                              %means temp2/dc picks the right proportion as
                              %well.
                  
           
            sopt(:,je,jjj,j2,j,TL)=(ones(da,1)-(temp2./dc)').*...
                                  (a+blimitage(TL)+...
               tlamb*(ylife(dy*dhip*(TL-1)+dy*(j2-1)+jjj)*ye(je,TL)+r*a).^(1-ttau)+beq)-...
               blimitage(TL);
           
         
        
      
                
         end % for j=1:dbet, j is for bet
      end % for je
      end % for jjj=1:dy, jjj is for h_t
     end %for j2=1:dhip, j2 is for hip group

      clear vri;

      %%%%%%%%%%%% working people, all other ages

      for j0=TL-1:-1:1 % j0 is the age index
          
          for j2=1:dhip %j2 is for the hip group
              
         for jjj=1:dy   % jjj is the index for ht
             for je=1:dye

            % first compute the utility matrix for each value of consumption
            % (rows)and today's assets (columns)
            c=linspace(0,1,dc+1)';
            c=c(2:dc+1)*(a'+blimitage(j0)+tlamb*(...
                ylife(dy*dhip*(j0-1)+dy*(j2-1)+jjj)*ye(je,j0)+r*a').^(1-ttau)+beq);     
            c(c<0)=NaN;
            u=c.^(1-sig)/(1-sig);

            % compute tomorrow's assets when not inheriting (=savings)
            % for each possible value of consumption and assets today
            s=-c+ones(dc,1)*(a'+blimitage(j0)+tlamb*(...
                ylife(dy*dhip*(j0-1)+dy*(j2-1)+jjj)*ye(je,j0)+r*a').^(1-ttau)+beq)-...
                blimitage(j0);
            st=reshape(s,1,da*dc);

            for j=1:dbet % j is the index for beta

               % we need to get a "big" extreme point for the value funtion
               % to avoid problem when we will interpolate. consider here
               % the value function for the guy that inherits
               
               vti=zeros(da+dagg,dye,dy);
               
               for i=1:dye
                   
               vti(:,i,:)=[vt(:,i,:,j2,j,j0+1);...
                   permute(intpcf*permute(vt(da-dpol+1:da,i,:,j2,j,j0+1),[1 3 2]),[1 3 2])];               
               end
               % note j2 is constant across lifetime
               
               %now vti (dar?,dye,dy)
               % we want to reduce the second dimension
               % idd assumption makes this very easy
               
               vti=repmat(yedtb(:,j0+1)',da+dagg,1,dy).*vti;
               vti=permute(sum(vti,2),[1 3 2]);
               
               vti=vti*Qylife(jjj,:,j0)';
               

               % let's consider the guy who has not inherited and interpolate
               % her value function
               vst=interp1q(ai,vti,st')';
               vst=reshape(vst,dc,da);
               % put toghether various pieces and compute the value function
               % for the guy who has already inherited (jj=dy+1)
               [temp1, temp2]=...
               max(u+sur(j0)*bet(j)*vst);
               vt(:,je,jjj,j2,j,j0)=temp1';
               
               vt(isnan(vt(:,je,jjj,j2,j,j0)),je,jjj,j2,j,j0)=-1e5;

               %polt(:,je,jjj,j2,j,j0)=temp2';
               copt(:,je,jjj,j2,j,j0)=(temp2./dc)'.*(a+blimitage(j0)+...
                                     tlamb*(...
                                     ylife(dy*dhip*(j0-1)+dy*(j2-1)+jjj)*ye(je,j0)+r*a).^(1-ttau)+beq);
               sopt(:,je,jjj,j2,j,j0)=(ones(da,1)-(temp2./dc)').*...
               (a+blimitage(j0)+tlamb*(...
               ylife(dy*dhip*(j0-1)+dy*(j2-1)+jjj)*ye(je,j0)+r*a).^(1-ttau)+beq)-...
               blimitage(j0);
               

            end % for j=1:dbet, j is for bet
             end % for je=1:dye
         end % for jjj=1:dy, jjj is for h_t
          end %for j2=1:dhip, j2 is for hip group
      end % for j0, j0 is for age
      disp('value function computed');

      %%%%%%%%%%%%% the value function computation is done
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      clear  c cc evat  s st temp1 temp2 vat vrr vst vti;

    %% 6.4. Invariant distributions
      %%%%%%%%%%%%% compute the invariant distribution 
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      t1=clock;
      % Since we have agents that die and are born at each period, the rows
      % of the transition matrix M do not sum up to one (people that die
      % disappear from the system and at each period a new cohort is born):
      % the people's distribution over the state variables evolves according
      % to: n_{t+1}' = n_t' M + q_{t+1}'
      % where q_{t+1} is the distribution of the people that are born at t+1

      % how is M structured? the more "external" state variable is beta:
      % the blocks for different beta's are in the main diagonal of the
      % matrix M (note that beta stays constant over the agent's lifetime).
      % Then the second variable is the agent's age.
      % For the first TL ages we distinguish by hip group and productivity level.
      % Finally we distinguish by assets
      % For ages from TL+1 to T the agent is retired, therefore for each age
      % we only distinguish by assets (productivity at birth is not relevant because
      % they do not inherit any more and current productivity is not relevant
      % because they do not work any more but just receive pensions).
      % transition matrix initialization
      
      %Now, depending on how many states we have, use normal indexation
      %(faster, takes up more memory) or sparse (slower, takes less memory)
      %{
      if (TL*dhip*dy*da+(T-TL)*da)*dbet<40000,
      M=zeros((TL*dhip*dy*da+(T-TL)*da)*dbet,(TL*dhip*dy*da+(T-TL)*da)*dbet);
      else
      M=sparse((TL*dhip*dy*da+(T-TL)*da)*dbet,(TL*dhip*dy*da+(T-TL)*da)*dbet);
      end
      %}

      %%%%%%%% consider retired people
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % To compute M we need first to reconduce the optimal asset choice
      % given by the policy functions on the same grid for the assets on which
      % we keep track of the value function (which is coarser to save on memory
      % space).
      % reconduce the ar indexes to the a indexes: note that this
      % transformation and the analogous ones that will follow crucially
      % depend on how the grids are formed. The point is the following:
      % I want to find the index to which a certain level of a corresponds
      %
      % Use that a=f(i), we can infer i=f-1(a)
      %
      % Example a = (0:sqrt(maxa)).^2
      % then sqrt(a) = (0:sqrt(maxa))
      % then sqrt(a(i)) is a point in that grid, that corresponds to
      % sqrt(a(i))= sqrt(maxa)*(i-1)/(da-1) 
      % since the first point is i=1 > 0 and the last point is i=da >1
      % More in general, we could say
      % sqrt(a(i)) = initialpoint + (sqrt(maxa)-initialpoint)*(i-1)/(da-1)
      %
      % We now want to get i out of there, just rearrange
      %
      % i = sqrt(a(i))/sqrt(maxa)*(da-1) + 1
      %
      
         indexr=zeros(size(solr));
         indexr(solr>=0)=floor((da-azero)*sqrt(solr(solr>=0))/sqrt(maxa)+1)+azero-1;
         if min(a)<0
         indexr(solr<0)=max(azero-ceil((azero-1)*sqrt(-solr(solr<0))/sqrt(-min(a))+1)+1,1);
         end
         
         indexr=(indexr-da+1).*(indexr<da)+da-1;
         restr=max(min((solr-a(indexr))./(a(indexr+1)-a(indexr)),1),0);

      
      % check the endpoint on both grids
      if ar(dar)>a(da)
         error('you have a problem with the asset grids');
      end



      % now find the distance from the smallest closer element to the element
      % in sopt that we are considering. 
      % This is because we will approximate the transition matrix and
      % the corresponding invariant distribution using the grid a and
      % giving weight to the two level of assets between which each optimal
      % choice is included. These weight will be the distances form those
      % two point (the total distance is obviously normalized to one since
      % we are talking about weights).
      
      
      


      %% take the optimal policy function for retired people (which lies on the
      % ar grid) and convert each of its points in the corresponding smaller,
      % closer value of assets on the grid a (indexr) and its distance from it
      % (restr)

      % find the index corresponding to the optimal savings using the same
      % mapping described above for the retired agents.
      
         indext=zeros(size(sopt));
         indext(sopt>=0)=floor((da-azero)*sqrt(sopt(sopt>=0))/sqrt(maxa)+1)+azero-1;
         if min(a)<0
         indext(sopt<0)=max(azero-ceil((azero-1)*sqrt(-sopt(sopt<0))/sqrt(-min(a))+1)+1,1);
         end
      
      
      indext=(indext-da+1).*(indext<da)+da-1;

      % now find the distance form the smallest closer element to the
      % element in sopt that we are considering. Note the we put on it
      % an upper bound of 1
      % in the case we would go over the maximum allowed grid for a. This is
      % because we will approximate the transition matrix and the corresponding
      % invariant distribution using the grid a and giving weight to the two
      % level of assets between which each optimal choice is included. These
      % weight will be the distances form those two point (the total distance
      % is obviously normalized to one since we are talking about weights).
      restt=max(min((sopt-a(indext))./(a(indext+1)-a(indext)),1),0);
       
      
      % We want to store indices
      % and numbers and then create the sparse matrix. Array M2 will store
      % rows, array M3 will store columns and array M4 will store values.
      
      %This makes sparse matrix indexation faster. Then Matlab just reads:
      %store value M4 in row M2, column M3.
      
      %Now fill them in. All matrices are
      %2*(TL*dhip*dy^2*da+(T-TL*da))*dbet. The 2 is there because for each
      %observation today we need to fill two cells: lower asset point and
      %higher asset point. So basically the first
      %dbet*(TL*dhip*dy^2*da+(T-TL*da) is the lower asset point on which we
      %put weight given a certain position today, and the second one is the
      %higher one.
      
      %The other particular thing is the dy^2. Why do we need it? Well,
      %from each dy we can transition to any other dy tomorrow. Therefore,
      %we need for each dy (1:dy)*da gridpoints to fill. We do them
      %contiguously, i.e., for each dy we have a row of contiguous
      %dy*da possible destinations that we fill. They are orthogonal, so we
      %solve it directly with the relevant Qylife.
                  
      newM2=zeros(2*dbet*((TL-1)*dhip*dy*dy*dye*dye*da+dhip*dye*dy*da+(T-TL-1)*da*dy*dye),1);
      newM3=newM2;
      newM4=newM2;
      
      fulld=(TL-1)*dhip*dy*dy*dye*dye*da+dhip*dy*dye*da+(T-TL-1)*da*dy*dye;
      fulldm=TL*dhip*dy*dye*da+(T-TL-1)*da*dy*dye;
      fullw=(TL-1)*dhip*dy*dy*dye*dye*da+dhip*dy*dye*da;
      fullr=dbet*((TL-1)*dhip*dy*dy*dye*dye*da+dhip*dy*dye*da+(T-TL-1)*da*dy*dye);
      %size newM2 is 2*fullr
      
      %Start with retired people: keep them in same j2
      for j=1:dbet
        for jj=1:TR-1
          for j2=1:dy %now they also vary by previous earnings
              for j3=1:dye % and by previous transitory shock
           for jjj=1:da
              newM2((j-1)*(fulld)+fullw+...
              (jj-1)*da*dy*dye+(j2-1)*da*dye+(j3-1)*da+jjj)=(j-1)*(fulldm)+TL*dhip*dy*da*dye+...
              (jj-1)*da*dy*dye+(j2-1)*da*dye+(j3-1)*da+jjj;
              newM2(fullr+...
              (j-1)*(fulld)+fullw+...
              (jj-1)*da*dy*dye+(j2-1)*da*dye+(j3-1)*da+jjj)=(j-1)*(fulldm)+TL*dhip*dy*da*dye+...
              (jj-1)*da*dy*dye+(j2-1)*da*dye+(j3-1)*da+jjj;
              newM3((j-1)*(fulld)+fullw+...
              (jj-1)*da*dy*dye+(j2-1)*da*dye+(j3-1)*da+jjj)=(j-1)*(fulldm)+TL*dhip*dy*da*dye+jj*da*dy*dye+...
              (j2-1)*da*dye+(j3-1)*da+indexr(jjj,j2,j,jj);
              newM3(fullr+...
              (j-1)*(fulld)+fullw+...
              (jj-1)*da*dy*dye+(j2-1)*da*dye+(j3-1)*da+jjj)=(j-1)*(fulldm)+TL*dhip*dy*da*dye+jj*da*dy*dye+...
              (j2-1)*da*dye+(j3-1)*da+indexr(jjj,j2,j,jj)+1;
              newM4((j-1)*(fulld)+fullw+...
              (jj-1)*da*dy*dye+(j2-1)*da*dye+(j3-1)*da+jjj)=sur(jj+TL)*(1-restr(jjj,j2,j,jj));
          %first asset point: gets 1-restr
              newM4(fullr+...
              (j-1)*(fulld)+fullw+...
              (jj-1)*da*dy*dye+(j2-1)*da*dye+(j3-1)*da+jjj)=sur(jj+TL)*restr(jjj,j2,j,jj);
          %second asset point, gets restr
           end % for jjj=1:da,
              end
          end
        end % for jj=1:TR-1,
      end % for j=1:dbet;
      
      %Now working age population:
      
      
      for j=1:dbet
         for jj=1:TL-1
             for j2=1:dhip
            for jjj=1:dy
                for je=1:dye % EXPLOIT iid-ness. This means your current je
                    % does not matter for your future je. I loop over each
                    % je but distribute them equally tomorrow. I do not
                    % need to loop over je tomorrow if I index well
               for j4=1:da
                  newM2((j-1)*(fulld)+(jj-1)*dye*dye*dy*dy*dhip*da+...
                      (j2-1)*dye*dye*dy*dy*da+(jjj-1)*dye*dye*da*dy+...
                      (je-1)*dye*da*dy+(0:dy*dye-1)*da+j4)=...
                      repmat((j-1)*(fulldm)+(jj-1)*dye*dy*dhip*da+...
                      (j2-1)*dye*dy*da+(jjj-1)*dye*da+(je-1)*da+j4,dy*dye,1);
                  %Note the (0:dy-1)*da: we are filling all of those. M2 is
                  %just a declaration, so in the RHS of M2 we find what the
                  %ACTUAL row of M this is going to be.
                  
                  %Saying it in another way, notice newM2, M3, M4 are
                  %LARGER than what M is going to be in order to store
                  %several points that we are going to need in order to
                  %write M. So in the left hand side, between parenthesis,
                  %we will find dy*dy, while in the RHS we only find dy.
                  
                  
                  newM3((j-1)*(fulld)+(jj-1)*dye*dye*dy*dy*dhip*da+...
                      (j2-1)*dye*dye*dy*dy*da+(jjj-1)*dye*dye*dy*da+...
                      (je-1)*dye*da*dy+(0:dy*dye-1)*da+j4)=(j-...
                      1)*(fulldm)+jj*dye*dhip*dy*da+...
                  (j2-1)*dye*dy*da+...
                  (0:dy*dye-1)*da+indext(j4,je,jjj,j2,j,jj);
                  %In the RHS of M3 we find the actual COLUMN of M this is
                  %going to be.
                  
                  newM4((j-1)*(fulld)+(jj-1)*dye*dye*dy*dy*dhip*da+...
                      (j2-1)*dye*dye*dy*dy*da+(jjj-1)*dye*dye*dy*da+...
                      (je-1)*dye*da*dy+(0:dy*dye-1)*da+j4)=...
                      sur(jj)*kron(Qylife(jjj,:,jj),yedtb(:,jj)')*...
                  (1-restt(j4,je,jjj,j2,j,jj));
                  %Here we fill the relevant dy values.
                  
                  %now upper point of asset holdings:
                  newM2(fullr+...
                      (j-1)*(fulld)+(jj-1)*dye*dye*dy*dy*dhip*da+...
                      (j2-1)*dye*dye*dy*dy*da+(jjj-1)*dye*dye*da*dy+...
                      (je-1)*dye*da*dy+(0:dy*dye-1)*da+j4)=...
                      (j-1)*(fulldm)+(jj-1)*dye*dhip*dy*da+...
                      (j2-1)*dye*dy*da+(jjj-1)*dye*da+(je-1)*da+j4;
                  newM3(fullr+...
                      (j-1)*(fulld)+(jj-1)*dye*dye*dy*dy*dhip*da+...
                      (j2-1)*dye*dye*dy*dy*da+(jjj-1)*dye*dye*da*dy+...
                      (je-1)*dye*da*dy+(0:dye*dy-1)*da+j4)=...
                      (j-1)*(fulldm)+jj*dye*dhip*dy*da+...
                      (j2-1)*dye*dy*da+...
                  (0:dy*dye-1)*da+indext(j4,je,jjj,j2,j,jj)+1;
                  newM4(fullr+...
                      (j-1)*(fulld)+(jj-1)*dye*dye*dy*dy*dhip*da+...
                      (j2-1)*dye*dye*dy*dy*da+(jjj-1)*dye*dye*da*dy+...
                      (je-1)*dye*dy*da+(0:dye*dy-1)*da+j4)=...
                      sur(jj)*kron(Qylife(jjj,:,jj),yedtb(:,jj)')*...
                  restt(j4,je,jjj,j2,j,jj);
               end % for j4=1:da,
                end % for je=1:dye
            end % for jjj=1:da,
             end %for j2=1:dhip
         end % for jj=1:TR-1,
      end % for j=1:dbet;
      clear jj;

      %%%%%%%% fill in the blocks for the working people at TL.
      % Now they need to retire towards their same jjj
      for j=1:dbet
        for j2=1:dhip
         for jjj=1:dy
           for je=1:dye
            for j4=1:da
               % this attributes the relevant mass
               % to the closest smaller value for tomorrow's assets
              
               
               newM2((j-1)*(fulld)+(TL-1)*dye*dye*dhip*dy*dy*da+...
                   (j2-1)*dye*dy*da+(jjj-1)*dye*da+(je-1)*da+j4)=(j-1)*(fulldm)+...
                   (TL-1)*dhip*dye*dy*da+(j2-1)*dye*dy*da+(jjj-1)*dye*da+...
                   (je-1)*da+j4;
                   
               newM3((j-1)*(fulld)+(TL-1)*dhip*dye*dye*dy*dy*da+...
                   (j2-1)*dye*dy*da+(jjj-1)*dye*da+(je-1)*da+j4)=(j-1)*(fulldm)+...
                   TL*dhip*dy*dye*da+...
                    (jjj-1)*da*dye+(je-1)*da+indext(j4,je,jjj,j2,j,TL);
                   
               newM4((j-1)*(fulld)+(TL-1)*dhip*dye*dye*dy*dy*da+...
                   (j2-1)*dye*dy*da+(jjj-1)*dye*da+(je-1)*da+j4)=...
                   sur(TL)*(1-restt(j4,je,jjj,j2,j,TL));
               
               newM2(fullr+...
                   (j-1)*(fulld)+(TL-1)*dhip*dye*dye*dy*dy*da+...
                   (j2-1)*dye*dy*da+(jjj-1)*dye*da+(je-1)*da+j4)=...
                   (j-1)*(fulldm)+(TL-1)*dye*dhip*dy*da+...
                   (j2-1)*dye*dy*da+(jjj-1)*dye*da+(je-1)*da+j4;
               
               newM3(fullr+...
                   (j-1)*(fulld)+(TL-1)*dhip*dye*dye*dy*dy*da+...
                   (j2-1)*dye*dy*da+(jjj-1)*dye*da+(je-1)*da+j4)=...
                   (j-1)*(fulldm)+TL*dye*dhip*dy*da+...
               (jjj-1)*da*dye+(je-1)*da+indext(j4,je,jjj,j2,j,TL)+1;
               
               newM4(fullr+...
                   (j-1)*(fulld)+(TL-1)*dhip*dye*dye*dy*dy*da+...
                   (j2-1)*dy*dye*da+(jjj-1)*dye*da+(je-1)*da+j4)=...
                   sur(TL)*restt(j4,je,jjj,j2,j,TL);
               
            end % for j4=1:da,
           end % for je
         end % for jjj=1:dy,
        end %for j2=1:dhip,
      end % for j=1:dbet;
      
      
      M=sparse(newM2,newM3,newM4,(TL*dhip*dy*dye*da+(T-TL)*da*dy*dye)*dbet,....
          (TL*dhip*dy*dye*da+(T-TL)*da*dy*dye)*dbet);
      %Now we just declare the M matrix.
      
      clear newM2 newM3 newM4
      
      % q is going to be the base for the invariant distribution
      %
      q=zeros(dbet*(TL*dhip*dy*dye*da+(T-TL)*da*dy*dye),1);
      % we need to initialise it with the distribution over the newborns
      
      % begin with those that have not inherited yet and start them at the
      % first asset level.

      tempvec=azero+kron((0:(dbet-1))'*(TL*dhip*dy*da+(T-TL)*da*dy*dye),ones(dye*dhip*dy,1))+...
              kron(ones(dbet,1),(0:(dhip*dy*dye-1))'*da); %removed +1 because
          % now it is azero
              % the first kron is for the beta block, the second one is for
              % coupling productivity and productivity at birth
      %q(tempvec)=kron(invherhip,yedtb(:,1));
      
      if initialwealth==0
          
      q(tempvec)=kron(invherhip,yedtb(:,1));
      
      elseif initialwealth==2
          load Inputs/iniwa
          % Need to reconstruct earnings deciles given the yedtb, invherhip
          % thing
          
          dtb25=kron(invherhip,yedtb(:,1));
          earn25=kron(ylife(1:dy),ones(dye,1)).*kron(ones(dy,1),ye(:,1));
          [~,I]=sort(earn25);
          cumdtb25=cumsum(dtb25(I));
          decile25=ceil(10*cumdtb25);
          newInd=zeros(size(I));
          
          unsorted = 1:length(earn25);
          newInd(I) = unsorted;
          
          rdecile25=decile25(newInd);
          
          for k=1:dy
              for k2=1:dye
                  q(tempvec((k-1)*dye+k2)+(1:da)-1)=invherhip(k).*yedtb(k2,1).*...
                      iniw(:,rdecile25((k-1)*dye+k2));
              end
          end
          
          
        elseif initialwealth==4
          load Inputs/iniwb
          % Need to reconstruct earnings deciles given the yedtb, invherhip
          % thing
          
          dtb25=kron(invherhip,yedtb(:,1));
          earn25=kron(ylife(1:dy),ones(dye,1)).*kron(ones(dy,1),ye(:,1));
          [~,I]=sort(earn25);
          cumdtb25=cumsum(dtb25(I));
          decile25=ceil(10*cumdtb25);
          newInd=zeros(size(I));
          
          unsorted = 1:length(earn25);
          newInd(I) = unsorted;
          
          rdecile25=decile25(newInd);
          
          for k=1:dy
              for k2=1:dye
                  q(tempvec((k-1)*dye+k2)+(1:da)-1)=invherhip(k).*yedtb(k2,1).*...
                      iniw(:,rdecile25((k-1)*dye+k2));
              end
          end
          
      end
     
        
      % compute the invariant distribution
      % note that we have population growth: n_{t+1}'=n_t' M+q_{t+1},
      % where q_{t+1} = q' gpop^(t+1); hence
      % n_{t+1}'=n_t' M+q gpop^(t+1); divide by gpop^(t+1)
      % n_{t+1}'/gpop^(t+1)=n_t'/gpop^t  M/gpop + q
      % redefine n_t'/gpop^t=p_t'; rewrite
      % p_{t+1}'=p_t' M/gpop +q  use this relationship to compute the
      % invariant distribution
      epsi=1;
      invm=q; % this is just an initialization
      %qiter=1;
      
      upi=full(sum(M,2));
      if max(upi)>1
          error('Problem in M')
      end
      

      while epsi>1e-08
         invm1=invm;
         invm=M'*invm/gpop+q; 
         epsi=max(abs(invm1-invm));
         %disp(qiter)
         %qiter=qiter+1;
         %disp(sum(invm))
      end
      q=q./sum(invm);
   
      
      clear q newM2 newM3 newM4 %liberate memory
      invm=invm./sum(invm);
      invm=full(invm);
      sutime2=etime(clock,t1)/60;
      fprintf('time elapsed %.6f \n',sutime2);
      disp('invariant distibution');
   
    %% 6.5. Some aggregates
    
   % sum all the people with different betas together
   invmt=reshape(invm,length(invm)/dbet,dbet);
   if dbet>1
      invmt=sum(invmt,2);
   end



   %%%%%%%%%%%% compute the number of retirees (sum over people who are
   %%%%%%%%%%%% 65 and older).
   numret=sum(invmt(TL*dhip*dy*dye*da+1:TL*dhip*dy*dye*da+TR*da*dy*dye));

   %%%%%%%%%%%% compute average capital (multiply each asset level for
   %%%%%%%%%%%% the number of people owning it, working or retired)
   avgk=sum(kron(ones(TL*dhip*dy*dye,1),a).*invmt(1:TL*dhip*dy*dye*da))+...
        sum(kron(ones(TR*dy*dye,1),a).*invmt(TL*dhip*dy*dye*da+1:TL*dhip*dy*dye*da+TR*da*dy*dye)); 
   
   yelife=reshape(kron(ones(1,dy),kron(ye,ones(da,1))),da*dye*dy*TL,1);

    
   %%%%%%%%%%%% compute average effective earnings
   aveinc=sum(((kron(ylife,ones(da*dye,1))).*yelife).*...
       invmt(1:TL*dhip*dy*dye*da));
         
   


