%%%---------------- generate_datas_oil_tristan.m----------------%%
%%% Program to generate the datas to two factor model based in the estimated parameters from
%%% Schartz- Smith oil future prices.
%%% 
%%%-------------------------------------------------
%%%
%%%% author: Karina Barbosa April 20
%%%%-------------------------------------------------
%%%  Program call function : "generator_noise.m" in order to design a
%%%  random vector with a desired covariance matrix.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
N=10000;%weeks
randn('state',1)% change this if you want different realizations
%%% parameter from
%%% Schartz- Smith oil future prices.
xo=0.5;
E0=log(20.2)-xo;
uee=1.15/100;    %sita 1
k=1.49;
sigx=28.6/100;
sige=14.5/100;
phoxe=0.3;
lambx=15.7/100;     %sita 1
ue=-1.25/100;     %sita 1
T=[0.0833    0.7500    1.0833    1.4167  2];

dt=1/48; % weekly Generate 
%dt=1;
nt=length(T);
%% %%%%%%% model %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x0=[xo;E0];
G=[exp(-k*dt) 0; 0 1]; %the dinamic of the model
c=[0;ue*dt];
% F=[];
% d=[];

%%%%%%%%%%%%% Covariance Matrices %%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%

%V=diag([0.042  0.003 0.0001 0.004 0.005]); 

V=0.005*eye(5);

W=[-expm1(-2*k*dt)*sigx^2/(2*k) (1-exp(-k*dt))*(phoxe*sigx*sige)/k;(1-exp(-k*dt))*(phoxe*sigx*sige)/k sige^2*dt ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%-----Vector of parameters-----%%%%%%%%%%%%%%%%%%%%%%%%
%Vp=[k sigx sige phoxe 0.042  0.003 0.0001 0.004 0.005 ue uee lambx E0];

%%
%%%%%%%%%%%%%%%%%%% construct output matrices %%%%%%%%%%%%%%%%%%%%%%
F=[];
d=[];
for i=1:nt
F=[F; exp(-k*T(i)) 1];
at=(uee+sige^2/2)*T(i)-expm1(-2*k*T(i))*sigx^2/(4*k)-expm1(-k*T(i))*(phoxe*sigx*sige-lambx)/k+expm1(-k*T(i))*lambx/k;
d=[d; at]; 
end
F=F';
xt=x0;
xt2=x0;
wt=[];
vt=[];
wt=generator_noise(W,N); %%% generator_noise(W,N) program to generate   wt with N terms given covariance matriz. It can be improved
vt=generator_noise(V,N); 


%% Simulation loop

Kyt=[]; %Vector:  output estimation 
Kxt=[];  %Vector:  state estimation 


%sdt=sqrt(dt);

for i=1:N %Generating the data 
xt=c+G*xt+wt(:,i);
yt=d+F'*xt+vt(:,i);
Kyt=[Kyt yt];
Kxt=[Kxt xt];
end

%save gera_oil T Kyt Kxt
Weeks=[1:1:N];
figure(1) 
plot(Weeks,exp(Kxt(2,:)),Weeks,exp(Kxt(1,:))) %equilibriom estate
grid on
legend('Long-term','Short-term')
title('generate state')
xlabel('Weeks')

figure(2)
plot(Weeks,exp(Kxt(2,:)+Kxt(1,:)),'.-',Weeks,exp(Kyt(1,:)),Weeks,exp(Kyt(2,:)),Weeks,exp(Kyt(3,:)),Weeks,exp(Kyt(4,:)),Weeks,exp(Kyt(5,:)))
grid on 
title('Output')
legend('Spot','1m','9m','13m','17m','24m') 
%Kyta=[(Kxt(2,:)+Kxt(1,:)); Kyt];