function  result = modelsimulation(T, number_of_runs,seed,varargin)
rng(seed,'twister');        % Set random seed
rndseed = rng;
rng(rndseed);

%% Parameters of the model
Tq=64; %number of days per quarter
h=1/Tq;%frequency 
sig=1.4; %risk aversion coefficient
c1=0.9; %stock price effect in IS equation   
rbar=0.01; %real ineterst rate in ss
betd=1/(1+h*rbar); %discount factor daily frequency
et=1; %inverse Frisch elascticity of labor supply
thet=0.55; %quarterly calvo price rigidity
thetd=1-h*(1-thet); %daily calvo price rigidity
gamd=(sig+et)*(1-thetd)*(1-betd*thetd)/thetd; 
kap=0.001; %stock price cost effect in phillips curve
alph=0.2; %extrapolation parameter in heuristics 
phi=10; %intensity of choice
for i = 1:length(varargin)
                                     if strcmpi(varargin{i},'phi')
                                     phi = varargin{i+1};
                                     end
 end
zet=0.5; %memory coefficient real
kap_s=0.04; %extrapolation for chartists and reversed speed for fundamentalists in financial sector
d=0.95; %memory financial sector
e=30; %intensity of choice financial sector
for i = 1:length(varargin)
                                     if strcmpi(varargin{i},'e')
                                     e = varargin{i+1};
                                     end
 end
a=1; %excess demand impact on stock price
k=0.2; %impact hh stock demand on excess demand
c2=0.5 ; %sensitivity to real stock price in hh stock demand
c3=1; %sensitivity to nomninal interest rate in hh stock demand
g=0.5; %bias fundamental stock price to output
for i = 1:length(varargin)
                                     if strcmpi(varargin{i},'g')
                                     g = varargin{i+1};
                                     end
 end
rho_s=0.7;  %daily stock shock persistence 
delta_x=0.5; %taylor rule output
for i = 1:length(varargin)
                                     if strcmpi(varargin{i},'delta_x')
                                     delta_x = varargin{i+1};
                                     end
 end
delta_pi= 1.5; %taylor rule inflation
for i = 1:length(varargin)
                                     if strcmpi(varargin{i},'delta_pi')
                                     delta_pi = varargin{i+1};
                                     end
end
delta_s=0; %taylor rule stock
 for i = 1:length(varargin)
                                     if strcmpi(varargin{i},'delta_s')
                                     delta_s = varargin{i+1};
                                     end
 end
sig_s=0.01; %stock price shock standard deviation. 1-2 percentage point per days on average
for i = 1:length(varargin)
                                     if strcmpi(varargin{i},'sig_s')
                                     sig_s = varargin{i+1};
                                     end
 end

%% matrices
A=[1 0 c1+h/sig*delta_s;-gamd h kap;0 0 1+a*k*c2+a*k*c3*delta_s];
B=[1 h*(1/sig-c1) c1; 0 betd*h 0; 0 0 0];

%% start values
x=zeros(T, number_of_runs); %output daily
pi=zeros(T, number_of_runs); %inflation daily
s=zeros(T, number_of_runs); %stock price
sf=zeros(T, number_of_runs); %fundamental stock price


E_s=zeros(T, number_of_runs); % expectation stock price
Et_s_c= zeros(T, number_of_runs); % expectation stock chartists
Et_s_F= zeros(T, number_of_runs); % expectation stock fundamentalists
Eq_tar_x=zeros(T, number_of_runs); % expectation output targetrs
Eq_tar_pi=zeros(T, number_of_runs); % expectation inflation targetrs
Eq_sta_x=zeros(T, number_of_runs); % expectation output static
Eq_sta_pi=zeros(T, number_of_runs); % expectation inflation static
Eq_ext_x=zeros(T, number_of_runs); % expectation output extrapolators
Eq_ext_pi=zeros(T, number_of_runs); % expectation inflation extrapolators
Eq_x=zeros(T, number_of_runs); % expectation output 
Eq_pi=zeros(T, number_of_runs); % expectation inflation 


eps_s=zeros(T, number_of_runs); % shock stock price

i=zeros(T, number_of_runs); %interest rate daily

A_tar_pi=zeros(T, number_of_runs); %attractivity targeters inflation
A_sta_pi=zeros(T, number_of_runs); %attractivity static inflation
A_tar_x=zeros(T, number_of_runs); %attractivity targeters output
A_sta_x=zeros(T, number_of_runs); %attractivity static output
A_ext_x=zeros(T, number_of_runs); %attractivity extrapolator output
A_ext_pi=zeros(T, number_of_runs); %attractivity extrapolator inflation

A_f_s=zeros(T, number_of_runs); %attractivity funademntalists stock price
A_c_s=zeros(T, number_of_runs); %attractivity charters stock price
A_nt_s=zeros(T, number_of_runs); %attractivity non-traders stock price

xq= zeros(T, number_of_runs); %output quarterly value
piq= zeros(T, number_of_runs); %inflation quarterly value
pq= zeros(T, number_of_runs); %price level quarterly

D=zeros(T, number_of_runs); % HH demand stock
D_f=zeros(T, number_of_runs); % fundamentalists demand stock
D_c=zeros(T, number_of_runs); % chartists demand stock


omega_x_tar=zeros (T,number_of_runs); %weights output targeters
omega_pi_tar=zeros(T,number_of_runs); %weights inflation targeters
omega_x_ext=zeros(T, number_of_runs); %weights output extrapolators
omega_pi_ext=zeros(T, number_of_runs); %weights inflation extrapolators
omega_x_sta=zeros(T, number_of_runs); %weights output static
omega_pi_sta=zeros(T, number_of_runs); %weights inflation static

omega_s_c=zeros(T, number_of_runs); %weights stock price charters
omega_s_f=zeros(T, number_of_runs); %weights stock price  fundamentalists
omega_s_nt=zeros(T, number_of_runs); %weights stock price  non-traders

xqq=zeros(T, number_of_runs); %output quarterly freq
piqq=zeros(T, number_of_runs); %inflation quarterly freq


%% special starting values
omega_x_tar(1:128,:)=1/3; 
omega_pi_tar(1:128,:)=1/3;
omega_x_sta(1:128,:)=1/3;%omega_pi_sta(1:128,:)=1/3;
omega_pi_ext(1:128,:)=1/3;
omega_x_ext(1:128,:)=1/3;
omega_s_c(1:128,:)=1/3;
omega_s_f(1:128,:)=1/3;
omega_s_nt(1:128,:)=1/3; 


%% shocks;
Randomz=randn(T, number_of_runs);
  
 %% Generate data
        for MC=1:number_of_runs
% start at first day of q3:
for t = (129:T)
  %from the first day of the third quarter until the last day in the last
  %quarter 

        
%% (((shocks)))
eps_s(t,MC)  = rho_s*eps_s(t-1,MC) +sig_s*Randomz(t,MC); 


    
        
 %% fundamental value
 sf(t,MC)=g*xq(t-1,MC);
 
 
 %% Attractivity
A_tar_pi(t,MC) = zet*A_tar_pi(t-1,MC) - (Eq_tar_pi (t-2,MC)-piq(t-1,MC))^2; 
A_sta_pi(t,MC) = zet*A_sta_pi(t-1,MC) - (Eq_sta_pi (t-2,MC)-piq(t-1,MC))^2; 
A_tar_x(t,MC) = zet*A_tar_x(t-1,MC) - (Eq_tar_x (t-2,MC)-xq(t-1,MC))^2; 
A_sta_x(t,MC) = zet*A_sta_x(t-1,MC) - (Eq_sta_x (t-2,MC)-xq(t-1,MC))^2; 
A_ext_x(t,MC) = zet*A_ext_x(t-1,MC) - (Eq_ext_x (t-2,MC)-xq(t-1,MC))^2; 
A_ext_pi(t,MC) = zet*A_ext_pi(t-1,MC) - (Eq_ext_pi (t-2,MC)-piq(t-1,MC))^2; 

A_f_s(t,MC) = d*A_f_s(t-1,MC) - (exp(s (t-2,MC))-exp(s(t-1,MC)))*D_f(t-2,MC); 
A_c_s(t,MC) = d*A_c_s(t-1,MC) - (exp(s (t-2,MC))-exp(s(t-1,MC)))*D_c(t-2,MC); 
A_nt_s(t,MC)=0;

%% (((weights)))
omega_x_tar(t,MC) = exp(phi*A_tar_x(t,MC))/(exp(phi * A_tar_x(t,MC)) + exp(phi * A_ext_x(t,MC))+ exp(phi * A_sta_x(t,MC)));
omega_pi_tar(t,MC)  = exp(phi*A_tar_pi(t,MC))/(exp(phi * A_tar_pi(t,MC)) + exp(phi * A_ext_pi(t,MC))+ exp(phi * A_sta_pi(t,MC)));
omega_x_ext(t,MC)  = exp(phi*A_ext_x(t,MC))/(exp(phi * A_tar_x(t,MC)) + exp(phi * A_ext_x(t,MC))+ exp(phi * A_sta_x(t,MC)));
omega_pi_ext(t,MC)=exp(phi*A_ext_pi(t,MC))/(exp(phi * A_tar_pi(t,MC)) + exp(phi * A_ext_pi(t,MC))+ exp(phi * A_sta_pi(t,MC)));
omega_x_sta(t,MC)=exp(phi*A_sta_x(t,MC))/(exp(phi * A_tar_x(t,MC)) + exp(phi * A_ext_x(t,MC))+ exp(phi * A_sta_x(t,MC)));
omega_pi_sta(t,MC)=exp(phi*A_sta_pi(t,MC))/(exp(phi * A_tar_pi(t,MC)) + exp(phi * A_ext_pi(t,MC))+ exp(phi * A_sta_pi(t,MC)));

omega_s_nt(t,MC)=exp(e*A_nt_s(t,MC))/(exp(e * A_f_s(t,MC)) + exp(e * A_c_s(t,MC))+ exp(e * A_nt_s(t,MC)));
omega_s_c(t,MC)=exp(e*A_c_s(t,MC))/(exp(e * A_f_s(t,MC)) + exp(e * A_c_s(t,MC))+ exp(e * A_nt_s(t,MC)));
omega_s_f(t,MC)=exp(e*A_f_s(t,MC))/(exp(e * A_f_s(t,MC)) + exp(e * A_c_s(t,MC))+ exp(e * A_nt_s(t,MC)));
%% Expectations

Eq_tar_x(t,MC)=0;
Eq_sta_x(t,MC)=xq(t-1,MC);
Eq_ext_x(t,MC)=xq(t-1,MC)+alph*(xq(t-1,MC)-xq(t-64,MC));
Eq_tar_pi(t,MC)=0;
Eq_sta_pi(t,MC)=piq(t-1,MC);
Eq_ext_pi(t,MC)=piq(t-1,MC)+alph*(piq(t-1,MC)-piq(t-64,MC));
Et_s_c(t,MC)=s(t-1,MC)+kap_s*(s(t-1,MC)-s(t-2,MC));
Et_s_F(t,MC)=s(t-1,MC)+kap_s*(sf(t-1,MC)-s(t-2,MC));
    Eq_x(t,MC) =   omega_x_tar(t,MC) *Eq_tar_x(t,MC)+omega_x_ext(t,MC) *Eq_ext_x(t,MC)+omega_x_sta(t,MC) *Eq_sta_x(t,MC);  % market expectation of output gap
    Eq_pi(t,MC) =  omega_pi_tar(t,MC) *Eq_tar_pi(t,MC)+omega_pi_ext(t,MC) *Eq_ext_pi(t,MC)+omega_pi_sta(t,MC) *Eq_sta_pi(t,MC);   % market expectation of inflation
    E_s(t,MC)=omega_s_c(t,MC) *Et_s_c(t,MC)+omega_s_f(t,MC) *Et_s_F(t,MC); 
 
    %% demand of stock
    D_c(t,MC)=(Et_s_c(t,MC)-s(t-1,MC));
    D_f(t,MC)=(Et_s_F(t,MC)-s(t-1,MC));
    D(t,MC)=xq(t,MC)-c2*(s(t-1,MC)-pq(t,MC))-c3*i(t-1,MC);
%% State Space Representation
Y = [Eq_x(t,MC); Eq_pi(t,MC); E_s(t,MC)];
Z = [-h/sig *piq(t,MC)*delta_pi-h/sig *(delta_x-g*delta_s)*xq(t,MC);0; ...
    s(t-1,MC)*(1+a*k*c2+a*k*c3*delta_s)+a*omega_s_c(t,MC)*D_c(t,MC)+omega_s_f(t,MC)*D_f(t,MC)+eps_s(t,MC)];
Model = B*Y + Z;
X = A\Model;
x(t,MC)= X(1);
pi(t,MC)= X(2);
s(t,MC)=X(3);

%% A way to calculate yq and piq correctly is to preserve the daily
    % timeline, making sure that the value only changes once per quarter..
    if mod(t,64) == 0  %last day of a quarter 
        xq(t,MC) = mean(x(t-63:t,MC) );
        piq(t,MC) = mean( pi(t-63:t,MC) ); 
        pq(t,MC)=sum(piq(64:64:t,MC));

    else
        xq(t,MC) = xq(t-1,MC);
        piq(t,MC)= piq(t-1,MC);
        pq(t,MC)=pq(t-1,MC);
    end
    

    
    %% Taylor rule 
i(t,MC)= delta_pi*piq(t,MC)+delta_x*xq(t,MC)+delta_s*s(t,MC); %daily interest rate change only once quarterly if delta_s=0

% Generate quarterly results (with quarterly values ONLY)
xqq(1:T/64,MC) = xq(64:64:T,MC);
piqq(1:T/64,MC) = piq(64:64:T,MC);

end
        end
%% %% Put all results into "result" variable
result.x  = x(129:T,1:number_of_runs);
result.Eq_sta_x  = Eq_sta_x(129:T,1:number_of_runs);


result.xq = xq(129:T,1:number_of_runs);
result.pi = pi(129:T,1:number_of_runs);
result.piq = piq(129:T,1:number_of_runs);
result.i  = i(129:T,1:number_of_runs);
result.Ex = Eq_x(129:T,1:number_of_runs);
result.Es = E_s(129:T,1:number_of_runs);
result.Esc = Et_s_c(129:T,1:number_of_runs);
result.Dc = D_c(129:T,1:number_of_runs);
result.Df = D_f(129:T,1:number_of_runs);

result.Dfomegaf =D_f(129:T,1:number_of_runs).*omega_s_f(129:T,1:number_of_runs);
result.Dcomegac =D_c(129:T,1:number_of_runs).*omega_s_c(129:T,1:number_of_runs);

result.Escomegac =  omega_s_c(129:T,1:number_of_runs).*Et_s_c(129:T,1:number_of_runs);
result.Esfomegaf =  omega_s_f(129:T,1:number_of_runs).*Et_s_F(129:T,1:number_of_runs);

result.Esf = Et_s_F(129:T,1:number_of_runs);
result.sf = sf(129:T,1:number_of_runs);
result.deltaEs = E_s(129:T,1:number_of_runs)-s(129:T,1:number_of_runs);
result.omega_x_tar = omega_x_tar(129:T,1:number_of_runs);
result.omega_x_sta = omega_x_sta(129:T,1:number_of_runs);
result.omega_x_ext = omega_x_ext(129:T,1:number_of_runs);
result.Epi = Eq_pi(129:T,1:number_of_runs);
result.omega_pi_tar = omega_pi_tar(129:T,1:number_of_runs);
result.omega_pi_sta = omega_pi_sta(129:T,1:number_of_runs);
result.omega_pi_ext = omega_pi_ext(129:T,1:number_of_runs);
result.A_pi_tar=A_tar_pi(129:T,1:number_of_runs);
result.A_pi_sta=A_sta_pi(129:T,1:number_of_runs);
result.A_x_sta=A_sta_x(129:T,1:number_of_runs);
result.A_x_tar=A_tar_x(129:T,1:number_of_runs);
result.A_s_c=A_c_s(129:T,1:number_of_runs);
result.A_s_f=A_f_s(129:T,1:number_of_runs);
result.s=s(129:T,1:number_of_runs);
result.omega_s_c=omega_s_c(129:T,1:number_of_runs);
result.omega_s_f=omega_s_f(129:T,1:number_of_runs);
result.omega_s_nt=omega_s_nt(129:T,1:number_of_runs);
result.epss=eps_s (129:T,1:number_of_runs);
result.xqq = xqq(3:T/64,1:number_of_runs);
result.piqq= piqq(3:T/64,1:number_of_runs);

result.varx=var(result.x);
result.varxqq=var(result.xqq);
result.varpi=var(result.pi);
result.varpiqq=var(result.piqq);
result.vars=var(result.s);
result.varr=var(result.i);

result.var_x=mean(result.varx);
result.var_s=mean(result.vars);
result.var_xqq=mean(result.varxqq);
result.var_pi=mean(result.varpi);
result.var_piqq=mean(result.varpiqq);


