clear
load plspec;
echo
%
%       This program is used to solve the following problem:
%    
%       (P)   max E sum{ beta^t*u(c1(t),..,cN(t),v2(t),param2) }
%          
%             s.t. c1(t)+i1(t)*IC(:,1) <= f1(v11(t),param11)
%                           :                  :
%                           :                  :
%                  cN(t)+iN(t)*IC(:,N) <= fN(v1N(t),param1N)
%                  X(t+1) = h(X(t),u(t))+g*epsilon(t)
echo
if info==2;
  disp('%                  y(t)=C*X(t)+w1(t)  (observables)')
  disp('%                                                                ')
end;
if info==3;
  disp('%                  u(t)=G*u1(t)+H*u2(t)     }              ')
  disp('%                  y1(t)=C1*X(t)+w1(t)      }  Information ')
  disp('%                  y2(t)=C2*X(t)+w2(t)      }    Structure')
  disp('%                                                                ')
end;
[rgetv,cgetv]=size(getv);
if (cgetv>10 & sum(getv(6,11:cgetv))>0);
  disp('The program, as written, assumes at most five production tech- ')
  disp('nologies (N=5) for printing results.  To get the results for cj,')
  disp('fj(.), j>5, that you have requested, you will have to edit    ')
  disp('"planner.m".  If the program continues, no error will result. ')
end;
if cgetv < 10;
  getv=[getv zeros(6,10-cgetv)];      
end;
param=[param1; param2];
rand('normal');
eval(desc); 
if (sol(2)==2 | sol(2)==3);
  disp('% ')
  disp('%    Below, we will require calculation of steady-state values of the')
  disp('%  nonstochastic version of our model.  These values are computed   ')
  disp('%  first to save time if poor starting values have been chosen.     ')
  disp('% ') 
end;
if ((nlcon==0 & sol(2)==3) | (nlcon==1 & sol(2)==2)); 
  error('There is an inconsistency in nlcon and sol')
end
if nlcon==0;
  if sol(2)==1;
    eval(ssxu);
  else;
    if overlap==0;
      [X,u,mmu]=steady(beta,A,B,D,E,param,crit,dampval,maxit,X0,u0,mu0);
    else;
      [X,u,mmu]=steadyst(beta,A,B,D,E,param,overlap,crit,dampval,maxit, ..
                         X0,u0,mu0);
    end;   
    x=D*X+E*u;
  end;
else;
  if sol(2)==1;
    eval(ssxu);
    adm=diag(.01*perdev1.*[X;u]);
    arg=find(abs([X;u])<.00001);
    for j=1:length(arg);
      adm(arg(j),arg(j))=.001;
    end;
    ns=length(X);
    nc=length(u);
    for j=1:ns+nc;
      jac(j,:)=(hfn([X;u]+adm(:,j),param3)'-hfn([X;u]-adm(:,j),param3)') ..
                /(2*adm(j,j));
    end;
    B=jac(ns+1:ns+nc,:)';
    A0=X-jac(1:ns,:)'*X-B*u;
    A=[1        zeros(1,ns-1);
       A0(2:ns)+jac(1,2:ns)' jac(2:ns,2:ns)'];
  else;
    [X,u,mmu,A,B]=nlsteady(beta,D,E,param,param3,crit,dampval,maxit,perdev1, ..
                         X0,u0,mu0);
    x=D*X+E*u;
  end;
end;
nx=length(x);
ns=length(X);
nc=length(u);
if sol(1)==1;
disp('%                                                                     ')
disp('%      Having specified the model and parameter values, we want to map')
disp('%  our (possibly) nonlinear problem into a linear-quadratic model which')
disp('%  is easily solved using the optimal linear regulator framework and, ')
disp('%  if states are hidden, the kalman filter.   That is, we want to map ')
disp('%  (P) into the following problem:                                    ')
disp('%                                                                     ')
disp('%  (P1) max E sum{ beta^t*(X(t)*Q*X(t)+u(t)*R*u(t)+2*X(t)*W*u(t) }    ')
disp('%                                                                     ')
disp('%          st.    X(t+1)=A*X(t)+B*u(t)+g*epsilon(t)                   ')
if info==2;
  disp('%                 y(t)=C*X(t)+w1(t)   y observed, X hidden ')
end;
if info==3;
  disp('%                 u(t)=G*u1(t)+H*u2(t)     u1,u2 chosen after seeing ')
  disp('%                 y1(t)=C1*X(t)+w1(t)      y1,y2, respectively ')
  disp('%                 y2(t)=C2*X(t)+w2(t)                               ')
end;
disp('%                                                                     ')
disp('%  where X is the state vector, u is the control vector.  The answer  ')
disp('%  to this problem is of the form:                                    ')
disp('%                                                                     ')
if info==1;
disp('%                       u(t)=-F*X(t).                                 ')
elseif info==2;
disp('%                u(t)=-F*Xhat(t),   Xhat(t)=E[X(t)|I(t)]              ')
else;
disp('%                     u1(t)=-F1*Xhat1(t),                             ')
disp('%                     u2(t)=-F2*Xhat1(t)-F3*u1(t),                    ')
disp('%  where                                                              ')
disp('%                     Xhat1(t)=E[X(t)|y1(t),I(t-1)]               ')
disp('%                     Xhat2(t)=E[X(t)|y2(t),y1(t),I(t-1)]       ')
disp('%                                                                     ')
end;
if (info==2 | info==3);
disp('%  Since our new objective function is linear-quadratic and our       ')
disp('%  constraints are linear, we can take advantage of the "separation   ')
disp('%  principle" which allows us to solve the prediction and control     ')
disp('%  problems independently.  The control problem is solved using the   ')
disp('%  optimal linear regulator framework, and the prediction is solved   ')
  if info==2; 
    disp('%  using a time invariant Kalman filter. ');
  else;
    disp('%  using a time-varying Kalman filter where the gain and variance ')
    disp('%  matrices take on two possible values depending on whether the  ')
    disp('%  time period is even or odd.                                    ')
  end;
end;
disp('%                                                                     ')
disp('%       The linear-quadratic utility function is found by taking a    ')
disp('%  second-order approximation of u(.) around the steady-state of the  ')
disp('%  nonstochastic version of the model.  With income less investment   ')
disp('%  substituted for consumption(s) in the utility function, the steady ')
disp('%  state vector, x, will include all of the production function       ')
disp('%  arguments and all of the utility function arguments with the       ')
disp('%  exception of c (i.e. [v11 .. v1N i1 .. iN v2], where there is no   ')
disp('%  duplication of variables ).  X and u will include steady state     ')
disp('%  values for states and controls, respectively.                      ')
disp('%                                                                     ')
disp('%       The stationary point of the model you chose is:               ')
disp('%                                                                     ')
x
disp(xstr);
X
disp(Xstr);
u
disp(ustr);
if sol(3)==1;
disp('%                                                                     ')
disp('%        You have chosen to approximate u(f1(.)-i1,...,v2,param2)     ')
disp('%  around the steady state using a second-order Taylor expansion.     ')
disp('%  Thus, it is necessary to evaluate the jacobian and hessian of u(.) ')
disp('%  at the stationary point.                                           ')
disp('%                                                                     ')
[ux,b,q]=ufn(x,param);
b=b(:);
q=q/2;           
else;
disp('%                                                                     ')
disp('%        You have chosen to approximate u(f1(.)-i1,...,v2,param2)     ')
disp('%  around the steady state using a method that approximates the       ')
disp('%  derivatives of u(.) without evaluating the first and second        ')
disp('%  derivatives.  Instead, the utility function must be evaluted at    ')
disp('%  points around the steady state.  For algebraic details, see        ')
disp('%  Kydland and Prescott (1982).                                       ')
disp('%                                                                     ')
adm=diag(.01*perdev2.*x);
arg=find(abs(x)<.00001);
for j=1:length(arg);
  adm(arg(j),arg(j))=.001;
end;
for j=1:nx
  b(j,1)=(ufn(x+adm(:,j),param)-ufn(x-adm(:,j),param))/(2*adm(j,j));
  q(j,j)=(ufn(x+adm(:,j),param)-2*ufn(x,param)+ufn(x-adm(:,j),param)) ...
         /(2*adm(j,j)^2);
  for k=1:j-1;
    q(j,k)=(ufn(x+adm(:,j)+adm(:,k),param)-ufn(x+adm(:,j)-adm(:,k),param)- ...
           ufn(x-adm(:,j)+adm(:,k),param)+ufn(x-adm(:,j)-adm(:,k),param)) ...
           /(8*adm(j,j)*adm(k,k));
    q(k,j)=q(j,k);
  end;
end;
end;
disp('%                                                                     ') 
disp('%       At this point, we have a linear-quadratic approximation of our') 
disp('%  return function, u(.).  In order to solve an optimal linear        ')
disp('%  regulator problem, it is necessary to distinguish the "states" and ') 
disp('%  the "controls".   Let X(t) and u(t) denote the vectors of state and')
disp('%  control variables at time t, respectively.  Then, we can write     ')
disp('%  x(t), which includes all of the variables from f1(.),..,fN(.), and ') 
disp('%  u(.) (except c) as:                                                ')
disp('%                                                                     ')
disp('%                     x(t)=D*X(t)+E*u(t).                         (1) ')
disp('%                                                                     ')
disp('%  Above, we found b and q such that:                                 ')
disp('%                                                                     ')
disp('%  u(f1(.)-i1(t)*IC(:,1),..)=u(x)+b*(x(t)-x)+(x(t)-x)*q*(x(t)-x)  (2) ')
disp('%                                                                     ')
disp('%  Therefore, in order to get the return function in terms of X(t) and') 
disp('%  u(t), we can use (1),(2), where D and E have been specified by you.')
disp('%                                                                     ')
disp('%       Given D, E, b, and q, it is easy to calculate Q, R, and W,    ')
disp('%  the matrices in the optimal linear regulator return function. Since')
disp('%  A and B,  which appear in the law of motion for X(t), have already')
disp('%  been specified by you, we have all of the information necessary to ')
disp('%  solve the control problem (P1).  Using the program "olrp.m", we    ')
disp('%  can obtain F, the feedback rule, and P, the solution to the Riccati')
if info~=3; disp('%  equation.  This is done below. '); end;
if info==3;
disp('%  equation.  Given P, we can then solve for F1, F2, and F3.  We maxi-')
disp('%  mize our value function with respect to u2 taking u1 as given. Then')
disp('%  we maximize the value function (with u2 replaced by u2*) with      ')
disp('%  respect to u1.  The value function depends directly on P.  This is ')
disp('%  done below.                                                        ')
end;
Q0=zeros(ns);
Q0(1,1)=ufn(x,param)-b'*x+x'*q*x;
d=(b'-2*x'*q)';
Q1=zeros(ns);
Q1(1,:)=d'*D;
Q1(:,1)=(d'*D)';
Q1=Q1/2;
Q1(1,1)=2*Q1(1,1);
Q=Q0+Q1+D'*q*D;
R=E'*q*E;
Q2=zeros(ns,nc);
Q2(1,:)=(d'*E)/2;
W=Q2+D'*q*E;
A
B
Q
R
W
[F,P]=olrp(beta,A,B,Q,R,W);
if (rank(R)<nc & check==1);
  disp('% ');
  disp('%  Note from check:');
  disp('%  Double cannot be run if R is singular');
  disp('% ');
end;
if (rank(R)==nc & check==1);
  [Ff,Pp]=dricc(beta,A,B,Q,R,W);
  disp('% ');
  disp('%  Note from check:')
  disp('%  The differences between gain matrices (F) and Riccati solutions (P)')
  disp('%  derived from double.m and dricc.m are:')
  disp(F-Ff)
  disp(P-Pp)
  if max(max(abs(F-Ff)))>1e-4; 
    F=Ff; 
    P=Pp; 
    disp('%  Since max(max(abs(double gain-dricc gain)))>1e-4, I have used the')
    disp('%  answers from dricc.m.')
    disp('% ');
  end;
end;
if info==1;
  FF=F;
  AA=A-B*F;
end;
if (info==2 | info==3);
  hidden=hidden(:)';
  nh=max(size(hidden));
  ny=length(V21);
  hmat=zeros(ns,nh);
  hidden=sort(hidden);
  for i=1:nh;
    hmat(hidden(i),i)=1;
  end;
%
% new section--preliminary
%
  if nh==1; sha=hmat'*A; sh=hmat'; else; sha=sum(hmat'*A); sh=sum(hmat'); end;
  tem1=find((abs(sha)>0)-sh>0);
  hidden=[hidden tem1];
  Ins=eye(ns);
  for i=1:length(tem1);
    C1=[C1;Ins(tem1(i),:)];
    C2=[C2;Ins(tem1(i),:)];
  end;
  V21=[V21                  zeros(ny,length(tem1));
       zeros(length(tem1)) .0000000001*eye(length(tem1))];
  V22=[V22                  zeros(ny,length(tem1));
       zeros(length(tem1)) .0000000001*eye(length(tem1))];
  nh=max(size(hidden));
  ny=length(V21);
  hmat=zeros(ns,nh);
  hidden=sort(hidden);
  for i=1:nh;
    hmat(hidden(i),i)=1;
  end;
  hmat2=hmat*hmat';
  Ih2=eye(ns)-hmat2;
  Ah=hmat'*A*hmat;
  Bh=hmat'*B;
  Sh=hmat'*g*Sigma*g*hmat;
end;
if info==2;
  Ch=C*hmat;
  gV32C=hmat'*g*(V3/V2)*Ch;
  [K,S]=kfilter(Ah,Ch,Sh,V2,hmat'*g*V3);
  FF=F*[Ih2 hmat];
  if ytknown==1;
    scv=S*(Ch'/(Ch*S*Ch'+V2));
    AA=[A-B*F*Ih2                                         -B*F*hmat;
      scv*Ch*(hmat'*A-gV32C*hmat')+gV32C*hmat'-Bh*F*Ih2   -Bh*F*hmat+ ..
                                                  (eye(nh)-scv*Ch)*(Ah-gV32C)];
    % Note:
    % For details: see Anderson and Moore, Optimal Filtering, eqs (5.5), (5.6)
  else;
    AA=[A-B*F*Ih2    -B*F*hmat;
        K*Ch*hmat'    Ah-K*Ch-Bh*F*hmat];
    % Note:
    % For details: see Sargent's filtering notes
    %
  end;
end;
if info==3;
  C1h=C1*hmat;
  C2h=C2*hmat;
  [ko,ke,so,se]=kf2stage(eye(nh),Ah,C1h,C2h,zeros(nh),Sh,V21,V22, ..
                       zeros(nh,ny),zeros(nh,ny));
  F3=inv(H'*(R+beta*B'*P*B)*H);
  F2=F3*H'*(beta*B'*P*A+W');
  F3=F3*H'*(R+beta*B'*P*B)*G;
  F1=((G-H*F3)'*(R+beta*B'*P*B)*(G-H*F3))\(G-H*F3)'*(beta*B'*P*(A-B*H*F2)+ ..
                                                     W'-R*H*F2);
  [row,nc1]=size(G);
  scv=se*(C2h'/(C2h*se*C2h'+V22));
  csc=scv*C2h;     
  icsc=eye(nh)-csc;
  akco=eye(nh)-ko*C1h;
  akce=Ah-ke*C2h;
  FF=[((G-H*F3)*F1+H*F2)*Ih2 zeros(nc,nh) (G-H*F3)*F1*hmat H*F2*hmat];
  AA=[A                                          zeros(ns,3*nh);
      ke*C2h*hmat'                               zeros(nh) akce      zeros(nh);
     (akco*ke*C2h+ko*C1h*Ah)*hmat'               zeros(nh) akco*akce zeros(nh);
     (icsc*(akco*ke*C2h+ko*C1h*Ah)+csc*Ah)*hmat' zeros(nh) icsc*akco*akce   ..
                                                                     zeros(nh)];
  AA=AA-[B;Bh;Bh;Bh]*FF;
  %
  % For details: see Ellen McGrattan's notes
end;
[RIC,CIC]=size(IC);
disp('%                                                                     ')
disp('%        We are now finished solving the control problem.  If some    ')
disp('%  or all of the states are hidden, then it is also necessary to solve')
disp('%  a filtering problem via the Kalman filter.  For all of the possible')
disp('%  information structures, however, we obtain a system of the form:   ')
disp('%                                                                     ')
disp('%                    W(t+1)=AA*W(t)+eta(t)                            ')
disp('%                     u(t)=-FF*W(t)                                   ')
disp('%                                                                     ')
disp('%  In this case,                                                      ')
AA
FF
if info==1;
  disp('%  where W(t)=X(t).                                                 ')
elseif info==2;
  disp('%  where W(t)=[X(t);Xhat1].                                         ')
else;
  disp('%  where W(t)=[X(t);Xhat0;Xhat1;Xhat2] and Xhat0 is used for compu- ')
  disp('%  tational purposes.                                               ')
end;
disp('% ')
disp('%       Given AA, FF, and E[eta*etaT], we can simulate time series for')
disp('%  the variables in the model by using "dlsim.m".  Below the requested')
disp('%  statistics for these time series are printed.                      ')
disp('%                                                                     ')
%
%   Compute means of variables to be printed
%
cz=sum(sum(getv));
mvec=ones(1,cz);
if morevars==1; temp=morev(ones(1,cz),1,param); cz=cz+max(size(temp)); end;
cc=1;
for i=1:5;
  for j=1:cgetv;
    if (getv(i,j) & i< 4 & detr(cc)); mvec(cc)=X(j); end;
    if (getv(i,j) & i==4 & detr(cc)); mvec(cc)=u(j); end;
    if (getv(i,j) & i==5 & detr(cc)); mvec(cc)=x(j); end;
    if getv(i,j)==1; cc=cc+1; end;
  end;
end;
lenx=1;
for i=1:N;
  if getv(6,2*i-1);
    if (i==1 & detr(cc)); mvec(cc)=ffn1(x(lenx:lenx+nov1(i)-1),param11); end;
    if (i==2 & detr(cc)); mvec(cc)=ffn2(x(lenx:lenx+nov1(i)-1),param12); end;
    if (i==3 & detr(cc)); mvec(cc)=ffn3(x(lenx:lenx+nov1(i)-1),param13); end;
    if (i==4 & detr(cc)); mvec(cc)=ffn4(x(lenx:lenx+nov1(i)-1),param14); end;
    if (i==5 & detr(cc)); mvec(cc)=ffn5(x(lenx:lenx+nov1(i)-1),param15); end;
    cc=cc+1;
  end;
  if (getv(6,2*i-1) & getv(6,2*i)); 
    if detr(cc);
      if IC==0;
        mvec(cc)=mvec(cc-1);
      else;
        mvec(cc)=mvec(cc-1)-x(sum(nov1)+(i-1)*RIC+1:sum(nov1)+i*RIC)'*IC(:,i);
      end;
    end;
    cc=cc+1;
  end;
  lenx=lenx+nov1(i);
end;
[rmv,cmv]=size(mvec);
if morevars==1; mvec=[mvec morev(mvec,1,param)]; end;
for i=length(detr):-1:cmv+1; 
  if detr(i)==0; mvec(i)=1; end;
end;
if min(abs(mvec))<1e-5; 
  error('Do not detrend variables with steady state values equal 0.');
end;
msd=zeros(1,cz);
mmean=zeros(1,cz);
sdsd=zeros(1,cz);
sdmean=zeros(1,cz);
mcor=zeros(N,cz);
sdcor=zeros(N,cz);
uservals=0;
gsg=g*Sigma*g';
if info==1;
  XI=X;
  vcov=gsg;
elseif info==2;
  XI=[X;hmat'*X];
  if ytknown==1;
    cross=gsg*hmat*Ch'*scv'+g*(V3/V2)*V3'*g'*hmat*(eye(nh)-scv*Ch)';
    vcov=[gsg           cross;
          zeros(nh,ns)  scv*Ch*Sh*Ch'*scv'+(eye(nh)-scv*Ch)*hmat'*g* ..
                        (V3/V2)*V3'*g'*hmat*(eye(nh)-scv*Ch)'+scv*V2*scv'+2* ..
                        scv*Ch*hmat'*g*(V3/V2)*V3'*g'*hmat*(eye(nh)-scv*Ch)'];
  else;
    vcov=[gsg             g*V3*K';
          zeros(nh,ns)    K*V2*K'];
  end;
else;      
  XI=[X;hmat'*X;hmat'*X;hmat'*X];
  zhs=zeros(nh,ns);
  zhh=zeros(nh);
  vcov=[gsg zhs'            gsg*hmat*C1h'*ko'    gsg*hmat*(csc'+C1h'*ko'*icsc');
        zhs ke*V22*ke'      ke*V22*ke'*akco'     ke*V22*ke'*akco'*icsc';
        zhs zhh  akco*ke*V22*ke'*akco'+ko*C1h*Sh*C1h'*ko'+ko*V21*ko' akco*ke* ..
           V22*ke'*akco'*icsc'+ko*C1h*Sh*(csc'+C1h'*ko'*icsc')+ko*V21*ko'*icsc';
        zhs zhh zhh icsc*akco*ke*V22*ke'*akco'*icsc'+(icsc*ko*C1h+csc)*Sh* ..
           (csc'+C1h'*ko'*icsc')+icsc*ko*V21*ko'*icsc'+scv*V22*scv'];
end;
vcov=triu(vcov)+triu(vcov,1)';
%
for i=1:nsim;
  clear z
  zt=zeros(lts,cz);
  lxi=length(XI);
  [ut,Xt]=dlsim(AA,vcov^.5,-FF,zeros(length(u),lxi),rand(lts,lxi),XI);
  xt=Xt(:,1:ns)*D'+ut*E';
  if info==1; Xhat1=Xt; end;
  if info==2; Xhat1=Xt(:,1:ns)*Ih2'+Xt(:,ns+1:ns+nh)*hmat'; end;
  if info==3; 
    Xhat1=Xt(:,1:ns)*Ih2'+Xt(:,ns+nh+1:ns+2*nh)*hmat';
    Xhat2=Xt(:,1:ns)*Ih2'+Xt(:,ns+2*nh+1:ns+3*nh)*hmat';
  else;
    Xhat2=Xt(:,1:ns);
  end;
  cc=1;
  for j=1:5;
    for k=1:cgetv;
      if (getv(j,k) & j==1); z(:,cc)=Xt(:,k); end;
      if (getv(j,k) & j==2); z(:,cc)=Xhat1(:,k); end;
      if (getv(j,k) & j==3); z(:,cc)=Xhat2(:,k); end;
      if (getv(j,k) & j==4); z(:,cc)=ut(:,k); end;
      if (getv(j,k) & j==5); z(:,cc)=xt(:,k); end;
      if getv(j,k)==1; cc=cc+1; end;
    end;
  end;
  clear Xt
  clear ut
  clear Xhat1
  clear Xhat2
  lenx=1;
  for j=1:N;
    if getv(6,2*j-1)==1;
      nof(2*j-1)=cc;
      if j==1; 
        for l=1:lts; z(l,cc)=ffn1(xt(l,lenx:lenx+nov1(1)-1)',param11); end;
      end;
      if j==2; 
        for l=1:lts; z(l,cc)=ffn2(xt(l,lenx:lenx+nov1(2)-1)',param12); end;
      end;
      if j==3;
        for l=1:lts; z(l,cc)=ffn3(xt(l,lenx:lenx+nov1(3)-1)',param13); end;
      end;
      if j==4;
        for l=1:lts; z(l,cc)=ffn4(xt(l,lenx:lenx+nov1(4)-1)',param14); end;
      end;
      if j==5;
        for l=1:lts; z(l,cc)=ffn5(xt(l,lenx:lenx+nov1(5)-1)',param15); end;
      end;
      cc=cc+1;
    end;
    if (getv(6,2*j-1) & getv(6,2*j));                  
      if IC==0; 
        z(:,cc)=z(:,cc-1);
      else;
        z(:,cc)=z(:,cc-1)-xt(:,sum(nov1)+(j-1)*RIC+1:sum(nov1)+j*RIC)*IC(:,j);
      end;
      cc=cc+1;
    end;
    lenx=lenx+nov1(j);
  end;
  if morevars==1;
    [zz zzstr]=morev(z,lts,param);
    z=[z zz];
  end;
  for j=1:cz; 
    if detr(j)==1; zt(:,j)=trend(z(:,j),lts,1600); end;
  end;
  dz=z-zt;
  if morefns==1; uservals=recurs(z,dz,uservals); end;
  msd=msd+(std(dz)./mvec)/nsim;
  mmean=mmean+mean(dz)/nsim;
  sdmean=sdmean+(mean(dz).^2)/nsim;
  sdsd=sdsd+((std(dz)./mvec).^2)/nsim;
  for k=1:N;
    if getv(6,2*k-1)==1; 
      ctemp=corr([dz(:,nof(2*k-1)) dz]); 
      mcor(k,:)=mcor(k,:)+ctemp(2:cz+1,1)'/nsim;
      sdcor(k,:)=sdcor(k,:)+((ctemp(2:cz+1,1)').^2)/nsim;
    end;
  end;
end;                % of i=1:nsim
sdsd=sqrt(sdsd-msd.^2);
sdmean=sqrt(sdmean-mmean.^2);
for k=1:N;
  sdcor(k,:)=sqrt(sdcor(k,:)-mcor(k,:).^2);
end;                
%
else;
X1=[X;u;mu];
disp('expath method is unfinished');
return
end;
%
disp('%                                                                     ')
disp('%  The means and standard deviations of the standard deviations of z *')
disp('%                                                                     ')
disp('%                                                                     ')
disp('  Means:   Standard Deviations: ')
disp([msd' sdsd']);
disp(zstr);
disp(zzstr);
disp(' ')
disp('* std(z-Az)/zbar, zbar=steady state value, Az=detrended z  (if detr=1) ')
disp('  std(z)   (if detr=0)                                                ')
disp(' ')
disp('%                                                                     ')
disp('%  The means and standard deviations of the correlations between      ')
disp('%  output(s) and elements of z  **                                    ')
for i=1:N;
disp('%                                                                     ')
disp('  Means:   Standard Deviations: ')
disp([mcor(i,:)' sdcor(i,:)']);
end;
disp(zstr);
disp(zzstr);
disp(' ')
disp('** corr([y-Ay z-Az]), Ay, Az=detrended y, z  (if detr=1)               ')
disp('   corr([y z])   (if detr=0)                                           ')
disp(' ')
disp(' ');
if morefns==1;
disp('  The following values were requested by the user via "RECURS.M":    ')
disp(uservals);
end;
%
disp('  To create plots, use matrix z saved in plan.mat.  I have not')
disp('  yet successfully automated plotting via batch.              ')
save plan
