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')