clear all;
close all;
clc;
T=40;  
lambda=0.4; beta=0.9;   teta=1.5;   alpha=1/2;   M1=1;   M2=2;
gamma=(lambda/(1-lambda))*(1-beta*(1-lambda))/(teta-(teta-1)*alpha);
delta=beta*(1-lambda);
Y=[ones(T/2,1)*(alpha*log(teta/(alpha*(teta-1)))+log(M1)); ones(T/2,1)*(alpha*log(teta/(alpha*(teta-1)))+log(M2))];
opt1=optimset('LargeScale','off','Display','off');
[P,fv]=fsolve(@sys,Y,opt1,Y(1),Y,T,gamma,delta);
lambda=0.6;
gamma=(lambda/(1-lambda))*(1-beta*(1-lambda))/(teta-(teta-1)*alpha);
delta=beta*(1-lambda);
Y2=[ones(T/2,1)*(alpha*log(teta/(alpha*(teta-1)))+log(M1)); ones(T/2,1)*(alpha*log(teta/(alpha*(teta-1)))+log(M2))];
[P2,fv]=fsolve(@sys,Y2,opt1,Y2(1),Y2,T,gamma,delta);
lambda=0.4; alpha=8/9;
gamma=(lambda/(1-lambda))*(1-beta*(1-lambda))/(teta-(teta-1)*alpha);
delta=beta*(1-lambda);
Y3=[ones(T/2,1)*(alpha*log(teta/(alpha*(teta-1)))+log(M1)); ones(T/2,1)*(alpha*log(teta/(alpha*(teta-1)))+log(M2))];
[P3,fv]=fsolve(@sys,Y3,opt1,Y3(1),Y3,T,gamma,delta);
alpha=1/2;  teta=1.1;
gamma=(lambda/(1-lambda))*(1-beta*(1-lambda))/(teta-(teta-1)*alpha);
delta=beta*(1-lambda);
Y4=[ones(T/2,1)*(alpha*log(teta/(alpha*(teta-1)))+log(M1)); ones(T/2,1)*(alpha*log(teta/(alpha*(teta-1)))+log(M2))];
[P4,fv]=fsolve(@sys,Y4,opt1,Y4(1),Y4,T,gamma,delta);
P2=P2.*Y./Y2;   P3=P3.*Y./Y3;   P4=P4.*Y./Y4;

t=1:T;
plot(t,exp(Y),t,exp(P),t,exp(P2),t,exp(P3),t,exp(P4),...
                'LineWidth',2)
legend('P*','P, baseline','P, \lambda > \lambdao','P, \alpha > \alphao','P, \theta < \thetao')
