% Prepare the price function of wheat without trend to plot it in R

load('../Results/WithoutTrend.mat')
a = Tableau{'Wheat','a'};
b = Tableau{'Wheat','b'};
delta = 0;
k = Tableau{'Wheat','k'};
r = 0.05;

invdemand           = @(d) a+b*d;

truncated_normal_rule(3,10,0,1,-5,5,fullfile('..','Results','GaussianQuadrature'))
model.shocks.e  = load(fullfile('..','Results','GaussianQuadrature_x.txt'));
model.shocks.w  = load(fullfile('..','Results','GaussianQuadrature_w.txt'));
TruncatedNormal = truncate(makedist('Normal'),-5,5);
model.funrand   = @(nrep,nper) TruncatedNormal.random(nrep,nper);
model.params    = [a b delta k r];

%% Define the grid of interpolation (equally space points)
interp.s           = linspace(-2,20,1000)';

options = struct('Display'     , 'iter',...
                 'InterpMethod', 'spline',...
                 'MaxIter'     , 1E3,...
                 'TolX'        , 1E-10);

%% Solve for the REE
interp       = SolveStorageDL(model,interp,options);

[Asim,Xsim,esim] = SimulStorage(model,interp,zeros(1E4,1),50+100);
vec = @(X) X(:);
Asim = vec(Asim(:,51:end));

s = interp.s;
x = interp.x;
Aobs = Aobs{:,'Wheat'};
Pobs = interp.cx{2}(Aobs);
pstar = Tableau{'Wheat','pstar'};
save(fullfile('..','Results','Wheat.mat'),'s','x','Aobs','Pobs','Asim','a','pstar','-v6')
