mostly field exp

Author: robertdkirkby

EvaluateFnOnAgentDist/FHorz/EvalFnOnAgentDist_AggVars_FHorz_Case1.m

@@ -84,6 +84,8 @@
     end
 end
 
+
+
 % If using e variable, do same for this
 if isfield(simoptions,'n_e')
     n_e=simoptions.n_e;
@@ -126,7 +128,7 @@
                     if all(size(simoptions.e_grid)==[sum(simoptions.n_e),1])
                         e_gridvals_J(:,:,jj)=gpuArray(CreateGridvals(simoptions.n_e,simoptions.e_grid,1));
                     else % already joint-grid
-                        e_gridvals_J(:,:,jj)=gpuArray(simoptions.e_grid,1);
+                        e_gridvals_J(:,:,jj)=gpuArray(simoptions.e_grid);
                     end
                 end
             else
@@ -135,7 +137,7 @@
                     if all(size(simoptions.e_grid)==[sum(simoptions.n_e),1])
                         e_gridvals_J(:,:,jj)=gpuArray(CreateGridvals(simoptions.n_e,simoptions.e_grid,1));
                     else % already joint-grid
-                        e_gridvals_J(:,:,jj)=gpuArray(simoptions.e_grid,1);
+                        e_gridvals_J(:,:,jj)=gpuArray(simoptions.e_grid);
                     end
                 end
             end
@@ -151,6 +153,8 @@
             n_z=[n_z,n_e];
             N_z=prod(n_z);
         end
+
+        simoptions=rmfield(simoptions,'n_e'); % Remove for subcommands as e is now in z
     end
 end
 

EvaluateFnOnAgentDist/FHorz/FieldExp/EvalFnOnAgentDist_AggVars_FHorz_Case1_FieldExp_PType.m

@@ -279,6 +279,7 @@
             end
             z_grid_J(:,jj)=gather(z_grid);
         end
+        simoptions_temp=rmfield(simoptions_temp,'ExogShockFn');
     else
         % This is just so that it makes things easier below because I can
         % take it as given that pi_z_J and z_grid_J exist
@@ -312,6 +313,7 @@
                 end
                 e_grid_J(:,jj)=gather(e_grid);
             end
+            simoptions_temp=rmfield(simoptions_temp,'EiidShockFn');
         else
             % This is just so that it makes things easier below because I can
             % take it as given that pi_e_J and e_grid_J exist
@@ -325,9 +327,10 @@
     simoptions_control_temp=simoptions_temp;
     simoptions_control_temp.z_grid_J=z_grid_J(:,TreatmentAgeRange(1):TreatmentAgeRange(2)+TreatmentDuration-1);
     if isfield(simoptions_temp,'n_e')
-        simoptions_control_temp.e_grid_J=e_grid_J(:,TreatmentAgeRange(1):TreatmentAgeRange(2)+TreatmentDuration-1);
+        simoptions_control_temp.e_grid=e_grid_J(:,TreatmentAgeRange(1):TreatmentAgeRange(2)+TreatmentDuration-1);
     end
 
+    
     %% Get the age weights out of control group as we need the later for treatment group, then reweight control group ages to reflect treatment duration
 
     % Find the age weights for the treatment group (get them from the marginal distribution of the control group over ages)
@@ -358,7 +361,6 @@
     StationaryDist_control_temp=reshape(StationaryDist_control_temp,[numel(StationaryDist_control_temp)/N_j_temp_FieldExp,N_j_temp_FieldExp]);
     StationaryDist_control_temp=StationaryDist_control_temp.*agereweight';
     StationaryDist_control_temp=StationaryDist_control_temp/sum(StationaryDist_control_temp(:));
-
 
     %% Calculate the aggregate variables for the control-group
     simoptions_control_temp.outputasstructure=0;

EvaluateFnOnAgentDist/FHorz/FieldExp/EvalFnOnAgentDist_AggVars_FHorz_Case1_FieldExp_Treatment.m

@@ -17,6 +17,7 @@
     simoptions.iterate=1;
     simoptions.tolerance=10^(-9);
     simoptions.outputkron=0; % If 1 then leave output in Kron form
+    simoptions.loopovere=0; % default is parallel over e, 1 will loop over e, 2 will parfor loop over e
 else
     %Check simoptions for missing fields, if there are some fill them with
     %the defaults
@@ -52,6 +53,9 @@
     if isfield(simoptions,'outputkron')==0
         simoptions.outputkron=0; % If 1 then leave output in Kron form
     end
+    if ~isfield(simoptions,'loopovere')
+        simoptions.loopovere=0; % default is parallel over e, 1 will loop over e, 2 will parfor loop over e
+    end
 end
 
 %%
@@ -82,7 +86,7 @@
 end
 
 if simoptions.parallel~=2 && simoptions.parallel~=4
-    Policy=gather(Policy);
+    PolicyKron=gather(PolicyKron);
     StationaryDist_Control=gather(StationaryDist_Control);    
     pi_z=gather(pi_z);
 end
@@ -158,10 +162,10 @@
 for j_p=TreatmentAgeRange(1):TreatmentAgeRange(2)
     % Pull the appropraite initial distribution of agents
     if N_ze==0
-        jequaloneDistKron=StationaryDist_Control(:,j_p);
+        jequaloneDistKron=reshape(StationaryDist_Control(:,j_p),[N_a,1]);
         PolicyKron_treat=PolicyKron(:,:,j_p:j_p+TreatmentDuration-1);
     else
-        jequaloneDistKron=StationaryDist_Control(:,:,j_p);
+        jequaloneDistKron=reshape(StationaryDist_Control(:,:,j_p),[N_a*N_ze,1]);
         if N_z==0 % just e
             PolicyKron_treat=PolicyKron(:,:,:,j_p:j_p+TreatmentDuration-1);
         elseif N_e==0 % Just z
@@ -192,11 +196,11 @@
             simoptions.pi_e_J=pi_e_J(:,j_p:j_p+TreatmentDuration);
         end
         if N_e==0
-            StationaryDist_treatment(:,:,:,j_p)=reshape(StationaryDist_FHorz_Case1_Iteration_raw(jequaloneDistKron,AgeWeightParamNames,PolicyKron_treat,N_d,N_a,N_z,TreatmentDuration,pi_z,Parameters,simoptions),[N_a,N_z,TreatmentDuration]);
+            StationaryDist_treatment(:,:,:,j_p)=reshape(StationaryDist_FHorz_Case1_Iteration_raw(jequaloneDistKron,AgeWeightParamNames,PolicyKron_treat,N_d,N_a,N_z,TreatmentDuration,pi_z_J,Parameters,simoptions),[N_a,N_z,TreatmentDuration]);
         elseif N_z==0
-            StationaryDist_treatment(:,:,:,j_p)=reshape(StationaryDist_FHorz_Case1_Iteration_noz_e_raw(jequaloneDistKron,AgeWeightParamNames,PolicyKron_treat,N_d,N_a,N_e,TreatmentDuration,simoptions.pi_e,Parameters,simoptions),[N_a,N_e,TreatmentDuration]); % e but no z
+            StationaryDist_treatment(:,:,:,j_p)=reshape(StationaryDist_FHorz_Case1_Iteration_noz_e_raw(jequaloneDistKron,AgeWeightParamNames,PolicyKron_treat,N_d,N_a,N_e,TreatmentDuration,pi_e_J,Parameters,simoptions),[N_a,N_e,TreatmentDuration]); % e but no z
         else
-            StationaryDist_treatment(:,:,:,j_p)=reshape(StationaryDist_FHorz_Case1_Iteration_e_raw(jequaloneDistKron,AgeWeightParamNames,PolicyKron_treat,N_d,N_a,N_z,N_e,TreatmentDuration,pi_z,simoptions.pi_e,Parameters,simoptions),[N_a,N_z*N_e,TreatmentDuration]);
+            StationaryDist_treatment(:,:,:,j_p)=reshape(StationaryDist_FHorz_Case1_Iteration_e_raw(jequaloneDistKron,AgeWeightParamNames,PolicyKron_treat,N_d,N_a,N_z,N_e,TreatmentDuration,pi_z_J,pi_e_J,Parameters,simoptions),[N_a,N_z*N_e,TreatmentDuration]);
         end
     end
 

StationaryDist/FHorz/FieldExp/StationaryDist_FHorz_Case1_FieldExp_Treatment.m

@@ -29,6 +29,7 @@
     simoptions.residualasset=0;
     simoptions.tolerance=10^(-9);
     simoptions.outputkron=0; % If 1 then leave output in Kron form
+    simoptions.loopovere=0; % default is parallel over e, 1 will loop over e, 2 will parfor loop over e
 else
     %Check simoptions for missing fields, if there are some fill them with
     %the defaults
@@ -80,6 +81,9 @@
     if ~isfield(simoptions,'outputkron')
         simoptions.outputkron=0; % If 1 then leave output in Kron form
     end
+    if ~isfield(simoptions,'loopovere')
+        simoptions.loopovere=0; % default is parallel over e, 1 will loop over e, 2 will parfor loop over e
+    end
 end
 
 

StationaryDist/FHorz/StationaryDist_FHorz_Case1.m

@@ -36,108 +36,6 @@
     StationaryDistKron=gpuArray(StationaryDistKron);
 end
 
-%% THIS COMMENTED OUT SECTION IS LEGACY CODE
-% Tried out doing this using cpu matrices, gpu matrices, sparse cpu
-% matrices, and sparce gpu matrices.
-% The sparse cpu matrices seem to be the fastest, and since they also use
-% least memory I have hardcoded them above.
-
-
-% if simoptions.parallel<2 || simoptions.parallel==3
-%     optaprime=gather(optaprime);
-%     jequaloneDistKron=gather(jequaloneDistKron);
-% end
-% 
-% if simoptions.parallel<2 % Full CPU
-% 
-%     StationaryDistKron=zeros(N_a,N_j);
-%     StationaryDistKron(:,1)=jequaloneDistKron;
-% 
-%     for jj=1:(N_j-1)
-% 
-%         if simoptions.verbose==1
-%             fprintf('Stationary Distribution iteration horizon: %i of %i \n',jj, N_j)
-%         end
-% 
-%         optaprime_jj=optaprime(1,:,jj);
-% 
-%         Gammatranspose=zeros(N_a,N_a); %probability of going to (a') given in (a)
-%         Gammatranspose(optaprime_jj+N_a*(0:1:N_a-1))=1;
-% 
-%         StationaryDistKron(:,jj+1)=Gammatranspose*StationaryDistKron(:,jj);
-%     end
-% 
-% elseif simoptions.parallel==2 % Full GPU
-% 
-%     StationaryDistKron=zeros(N_a,N_j,'gpuArray');
-%     StationaryDistKron(:,1)=jequaloneDistKron;
-% 
-%     for jj=1:(N_j-1)
-% 
-%         if simoptions.verbose==1
-%             fprintf('Stationary Distribution iteration horizon: %i of %i \n',jj, N_j)
-%         end
-% 
-%         optaprime_jj=optaprime(1,:,jj);
-% 
-%         Gammatranspose=zeros(N_a,N_a,'gpuArray');
-%         Gammatranspose(optaprime_jj+N_a*(gpuArray(0:1:N_a-1)))=1;
-% 
-%         StationaryDistKron(:,jj+1)=Gammatranspose*StationaryDistKron(:,jj);
-%     end
-% 
-% elseif simoptions.parallel==3 % Sparse CPU
-% 
-%     StationaryDistKron=zeros(N_a,N_j);
-%     StationaryDistKron(:,1)=jequaloneDistKron;
-% 
-%     StationaryDistKron_jj=sparse(jequaloneDistKron);
-% 
-%     for jj=1:(N_j-1)
-% 
-%         if simoptions.verbose==1
-%             fprintf('Stationary Distribution iteration horizon: %i of %i \n',jj, N_j)
-%         end
-% 
-%         optaprime_jj=optaprime(1,:,jj);
-% 
-%         Gammatranspose=sparse(N_a,N_a);
-%         Gammatranspose(optaprime_jj+N_a*(0:1:N_a-1))=1;
-% 
-%         StationaryDistKron_jj=Gammatranspose*StationaryDistKron_jj;
-%         StationaryDistKron(:,jj+1)=full(StationaryDistKron_jj);
-%     end
-% 
-%     if gpuDeviceCount>0 % Move the solution to the gpu if there is one
-%         StationaryDistKron=gpuArray(StationaryDistKron);
-%     end
-% 
-% elseif simoptions.parallel==4 % Sparse matrix instead of a standard matrix for P, on gpu
-% 
-%     StationaryDistKron=zeros(N_a,N_j,'gpuArray');
-%     StationaryDistKron(:,1)=jequaloneDistKron;
-% 
-%     StationaryDistKron_jj=sparse(StationaryDistKron(:,1));
-%     for jj=1:(N_j-1)
-% 
-%         if simoptions.verbose==1
-%             fprintf('Stationary Distribution iteration horizon: %i of %i \n',jj, N_j)
-%         end
-% 
-%         optaprime_jj=optaprime(1,:,jj);
-% 
-%         Gammatranspose=sparse(N_a,N_a);
-%         Gammatranspose(optaprime_jj+N_a*(0:1:N_a-1))=1;
-% 
-%         Gammatranspose=gpuArray(Gammatranspose);
-% 
-% %         StationaryDistKron(:,jj+1)=Gammatranspose*StationaryDistKron(:,jj); % Cannot index sparse gpuArray, so have to use StationaryDistKron_jj instead
-%         StationaryDistKron_jj=Gammatranspose*StationaryDistKron_jj;
-%         StationaryDistKron(:,jj+1)=full(StationaryDistKron_jj);
-%     end
-% 
-% end
-
 %%
 
 % Reweight the different ages based on 'AgeWeightParamNames'. (it is assumed there is only one Age Weight Parameter (name))