function [ys,check] = finfric158_steadystate(ys,junk);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The model code is based on Lawrence J. Christiano, Mathias Trabandt and Karl Walentin, 2011. Introducing financial frictions and unemployment into a small open economy model, Journal of Economic Dynamics and Control, Elsevier, vol. 35(12), pages 1999-2041.
% Adapted by Ginters Buss for Latvia as a small open economy in a monetary
% union.
% 
% If these codes are useful, cite: 
% Ginters Buss, 2016. "Financial frictions in Latvia," Empirical Economics, Springer, vol. 51(2), pages 547-575, September.
% 
% Copyright of this version: Latvijas Banka
% These codes are distributed by Latvijas Banka for academic purposes. 
% If in doubt, contact Latvijas Banka.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global M_ check;
persistent LastSuccessFsolve;

% Here we load the values of the deep parameters in a loop.
myNp = M_.param_nbr;                                            
for i = 1:myNp
    paramname = deblank(M_.param_names(i,:));
    eval([ paramname ' = M_.params(' int2str(i) ');']);
end

% set all working capital parameters to the one common estimated value
nuf=work_cap_para;
nustar=work_cap_para;
nux=work_cap_para;

%calculate no-indexation parameters
varkappad=1-kappad;
varkappax=1-kappax;
varkappamx=1-kappamx;
varkappami=1-kappami;
varkappamc=1-kappamc;
varkappaw=1-kappaw;

muz=1+muzPercent/100;
% get mupsi residually, so that muzplus is always the same (and must be given in .mod-file
mupsi=(muzplus/muz)^((1-alphaa)/alphaa);

weshare=wesharePercent/100;
 
sigmaL=sigmaLScaled*10;
Spp=SppScaled*10;   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Steady State Calculations%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Inflation rates
pic=pitarget;
pibar=pic;
pimx=pibar;
pimc=pibar;
pimi=pibar;
pii=pibar/mupsi;
piw=muzplus*pibar;
pix=pistar;


%Exchange rate growth and return on capital
s=pic/pistar;
%Rk=pibar*muzplus/betaa; %not valid in financial frictions model

% Nominal Interest Rate
R=(muzplus/betaa-taub)/(1-taub)*pibar;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%World Nominal Interest Rate
Rstar=(muzplus/betaa-taub)/(1-taub)*pistar;

%Working Capital Interest Rates
Rnustar=nustar*Rstar+1-nustar;
Rf=nuf*R+1-nuf;
Rx=nux*R+1-nux;

%Nonlinear pricing domestic intermediate
pitilded=pibar^kappad*pitarget^(1-kappad-varkappad)*pibreve^varkappad;

mc=1/lambdad*(1-betaa*xid*(pitilded/pibar)^(lambdad/(1-lambdad)))...
    /(1-betaa*xid*(pitilded/pibar)^(1/(1-lambdad)))...
    *((1-xid*(pitilded/pibar)^(1/(1-lambdad)))/(1-xid))^(1-lambdad);

phalo=((1-xid*(pitilded/pibar)^(1/(1-lambdad)))/(1-xid))^(1-lambdad)...
    /((1-xid*(pitilded/pibar)^(lambdad/(1-lambdad)))/(1-xid))^((1-lambdad)/lambdad);


%Nonlinear pricing exporters
pitildex=pix^kappax*pix^(1-kappax-varkappax)*pibreve^varkappax;

mcx=1/lambdax*((1-xix*(pitildex/pix)^(1/(1-lambdax)))/(1-xix))^((1-lambdax))...
    *(1-betaa*xix*(pitildex/pix)^(lambdax/(1-lambdax)))...
    /(1-betaa*xix*(pitildex/pix)^(1/(1-lambdax)));


phalox=((1-xix*(pitildex/pix)^(1/(1-lambdax)))/(1-xix))^(1-lambdax)...
    /((1-xix*(pitildex/pix)^(lambdax/(1-lambdax)))/(1-xix))^((1-lambdax)/lambdax);


%Nonlinear pricing consumption importer
pitildemc=pimc^kappamc*pitarget^(1-kappamc-varkappamc)*pibreve^varkappamc;

mcmc=1/lambdamc*((1-ximc*(pitildemc/pimc)^(1/(1-lambdamc)))/(1-ximc))^((1-lambdamc))...
    *(1-betaa*ximc*(pitildemc/pimc)^(lambdamc/(1-lambdamc)))...
    /(1-betaa*ximc*(pitildemc/pimc)^(1/(1-lambdamc)));

phalomc=((1-ximc*(pitildemc/pimc)^(1/(1-lambdamc)))/(1-ximc))^(1-lambdamc)...
    /((1-ximc*(pitildemc/pimc)^(lambdamc/(1-lambdamc)))/(1-ximc))^((1-lambdamc)/lambdamc);

%Nonlinear pricing investment importer
pitildemi=pimi^kappami*pitarget^(1-kappami-varkappami)*pibreve^varkappami;

mcmi=1/lambdami*((1-ximi*(pitildemi/pimi)^(1/(1-lambdami)))/(1-ximi))^((1-lambdami))...
    *(1-betaa*ximi*(pitildemi/pimi)^(lambdami/(1-lambdami)))...
    /(1-betaa*ximi*(pitildemi/pimi)^(1/(1-lambdami)));

phalomi=((1-ximi*(pitildemi/pimi)^(1/(1-lambdami)))/(1-ximi))^(1-lambdami)...
    /((1-ximi*(pitildemi/pimi)^(lambdami/(1-lambdami)))/(1-ximi))^((1-lambdami)/lambdami);

%Nonlinear pricing export importer
pitildemx=pimx^kappamx*pitarget^(1-kappamx-varkappamx)*pibreve^varkappamx;

mcmx=1/lambdamx*((1-ximx*(pitildemx/pimx)^(1/(1-lambdamx)))/(1-ximx))^((1-lambdamx))...
    *(1-betaa*ximx*(pitildemx/pimx)^(lambdamx/(1-lambdamx)))...
    /(1-betaa*ximx*(pitildemx/pimx)^(1/(1-lambdamx)));

phalomx=((1-ximx*(pitildemx/pimx)^(1/(1-lambdamx)))/(1-ximx))^(1-lambdamx)...
    /((1-ximx*(pitildemx/pimx)^(lambdamx/(1-lambdamx)))/(1-ximx))^((1-lambdamx)/lambdamx);

%Wage setting distortions and indexation
pitildew=pic^kappaw*pitarget^(1-kappaw-varkappaw)*pibreve^varkappaw*muzplus^thetaw;


whalo=((1-xiw*(pitildew/piw)^(1/(1-lambdaw)))/(1-xiw))^(1-lambdaw)...
    /((1-xiw*(pitildew/piw)^(lambdaw/(1-lambdaw)))/(1-xiw))^((1-lambdaw)/lambdaw);



%initial guess

% at mh mode 2010-02-15:
xxx=[1
     1
     1
    1  
    1
    1
    1
    1
    1
    %0.0270 %delta
    1];
    
max_fsolve_tries=10;
fsolve_tries=1;
check=1;

%use Fsolve to solve the core
while fsolve_tries<=max_fsolve_tries && check==1,

    if fsolve_tries == 1,
        if ~isempty(LastSuccessFsolve)
            xxx=LastSuccessFsolve;
        end

    elseif fsolve_tries>2
        ten=10;
        lowerbound=0.01;
   
        % below guesses from ctw_steadystate
        xxx(1)=1*(ten*rand*rand/2+lowerbound);%0.4580; %sigma
        xxx(2)=1*(ten*rand*rand/2+lowerbound);%0.3805; %omegabar
        xxx(3)=1*(ten*rand*rand/2+lowerbound);%1.0241;%Rk
        xxx(4)=1*(ten*rand*rand/2+lowerbound);%7.4321;%kbar
        xxx(5)=1*(ten*rand*rand/2+lowerbound);%1.1845;%wbar
        xxx(6)=1*(ten*rand*rand/2+lowerbound);%0.6773;%H
        xxx(7)=1*(ten*rand*rand/2+lowerbound);%0.1620;%ystar
        xxx(8)=1*(ten*rand*rand/2+lowerbound);%0.0471;%rkbar
        xxx(9)=1*(ten*rand*rand/2+lowerbound);%AL
        %xxx(10)=0.02*(ten*rand*rand/2+lowerbound);%delta
        xxx(10)=1*(ten*rand*rand/2+lowerbound);%varphi
        end      
     
   try

        opts=optimset('Display','off','MaxIter',200,'TolFun',1e-9,'TolX',1e-9);
        [yyy,fval,exitflag]=fsolve(@finfric_steadystatesolvecore_noDelta3,xxx,opts,pibar,mupsi,tauk,pkprime,muzplus,R,alphaa,Rf,epsilon,lambdad,...
            etag,Rx,etax,omegax,pmx,px,etaf,omegai,pinvest,etai,Upsilon,zetac,betaa,b,pc,tauc,zetah,sigmaL,omegac,...
            etac,Rnustar,tauw,tauy,etaa,taud,s,phalo,pitildew,lambdamc,pmc,mc,piw,phalomi,phalomc,...
            phalomx,Rstar,lambdami,pmi,lambdamx,phalox,lambdax,Rk,mu,we,Fomegabar,xiw,lambdaw,...
            intensity_target,iy_target,xy_target,nk_target,taumx,taumi,taumc,mcmx,mcmc,mcmi,taux,mcx,delta);
   

    fsolve_tries=fsolve_tries+1;

    if  exitflag~=1,
        %disp('Cannot find a steady state');
        ys=ones(size(ys))*1e+2;
        check=1;


    elseif exitflag==1
        sigma=yyy(1)*.41;
        omegabar=yyy(2)*.39;
        Rk=yyy(3)*1.05;
        kbar=yyy(4)*1.12;
        wbar=yyy(5)*.8;
        H=yyy(6)*.27;
        ystar=yyy(7)*.07;
        rkbar=yyy(8)*.13;
        AL=yyy(9)*4.16;
        %delta=yyy(10);
        varphi=yyy(10)*2.37;
        
        LastSuccessFsolve=yyy;
        
        %Relative import prices
pmx=taumx*Rnustar*varphi/mcmx;
pmc=taumc*Rnustar*varphi/mcmc;
pmi=taumi*Rnustar*varphi/mcmi;

%Relative consumption price
pc=(1-omegac+omegac*pmc^(1-etac))^(1/(1-etac));

%Relative investment price
pinvest=(1-omegai+omegai*pmi^(1-etai))^(1/(1-etai));

%Relative export price
px=(taux*Rx/(mcx*varphi))*(omegax*pmx^(1-etax)+1-omegax)^(1/(1-etax));

%Relative price of capital
pkprime=pinvest/Upsilon;


        %some abbreviations
        z       =   (log(omegabar)+sigma^2/2)/sigma;
        Gamma   =   normcdf(z-sigma)+omegabar*(1-normcdf(z));
        Gam_muG =   (1-mu)*normcdf(z-sigma)+omegabar*(1-normcdf(z));
        
        %ybar
        ybar=(phalo^(lambdad/(lambdad-1))*epsilon*(kbar/H/(mupsi*muzplus))^alphaa)...
            /(1+phalo^(lambdad/(lambdad-1))*(1/mc-phalo^(lambdad/(1-lambdad))))*H;
        
    
        %Further steady states
        nk=(pkprime*(1-Rk/R*Gam_muG));
        we=weshare*ybar;
        
        %use law of motion for net worth to calc gammaa
        gammaa=(1-we/nk/kbar)/((Rk-R-(Gamma-Gam_muG)*Rk)*pkprime/nk+R)*(pibar*muzplus);         
        n=nk*kbar;
        spread=omegabar*Rk/(1-n/(pkprime*kbar))-R;

        g=etag*ybar;
        q=varphi/pc;
        i=kbar*(1-(1-delta)/(muzplus*mupsi))/Upsilon;
        x=px^(-etaf)*ystar;
     
        c=1/(zetah*AL)*((1-xiw*(pitildew/piw)^(1/(1-lambdaw)))/(1-xiw))^(1-lambdaw*(1+sigmaL))...
            *(wbar*(1-tauy)/(lambdaw*(1+tauw)))*(1-betaa*xiw*(pitildew/piw)^(lambdaw*(1+sigmaL)/(1-lambdaw)))...
            *(zetac/(muzplus-b)*(muzplus-betaa*b)/(pc*(1+tauc)))/(1-betaa*xiw*(pitildew/piw)^(1/(1-lambdaw)))*H^(-sigmaL);
        %c=(zetac*(muzplus-betaa*b)/(pc*(1+tauc)*(muzplus-b))...
        %    /(zetah*AL*(H/(1-unemp))^sigmaL/wbar/((1-tauy)/(1+tauw))));
        phi=ybar*(1/mc-phalo^(lambdad/(1-lambdad)));
        a=etaa*ybar;
        k=kbar*u;
        psizplus=(zetac*(muzplus-betaa*b)/(muzplus-b))/(c*pc*(1+tauc));
        
        im=omegai*(pinvest/pmi)^(etai)*i;
        xm=omegax*(1/pmx*(omegax*pmx^(1-etax)+1-omegax)^(1/(1-etax)))^etax...
            *phalox^(-lambdax/(lambdax-1))*px^(-etaf)*ystar;
        cm=omegac*(pc/pmc)^(etac)*c;


        %nonlinear pricing variables
        Kd=lambdad*psizplus*ybar*mc/(1-betaa*xid*(pitilded/pibar)^(lambdad/(1-lambdad)));
        Fd=psizplus*ybar/(1-betaa*xid*(pitilded/pibar)^(1/(1-lambdad)));

        Kx=lambdax*psizplus*varphi*px*x*mcx/(1-betaa*xix*(pitildex/pix)^(lambdax/(1-lambdax)));
        Fx=psizplus*varphi*px*x/(1-betaa*xix*(pitildex/pix)^(1/(1-lambdax)));

        Kmc=lambdamc*psizplus*pmc*mcmc*cm/(1-betaa*ximc*(pitildemc/pimc)^(lambdamc/(1-lambdamc)));
        Fmc=psizplus*pmc*cm/(1-betaa*ximc*(pitildemc/pimc)^(1/(1-lambdamc)));

        Kmi=lambdami*psizplus*pmi*mcmi*im/(1-betaa*ximi*(pitildemi/pimi)^(lambdami/(1-lambdami)));
        Fmi=psizplus*pmi*im/(1-betaa*ximi*(pitildemi/pimi)^(1/(1-lambdami)));

        Kmx=lambdamx*psizplus*pmx*mcmx*xm/(1-betaa*ximx*(pitildemx/pimx)^(lambdamx/(1-lambdamx)));
        Fmx=psizplus*pmx*xm/(1-betaa*ximx*(pitildemx/pimx)^(1/(1-lambdamx)));

        
        Kw=(zetah*H^(1+sigmaL))/(1-betaa*xiw*(pitildew/piw)^(lambdaw*(1+sigmaL)/(1-lambdaw)));
        Fw=(psizplus/lambdaw*H*(1-tauy)/(1+tauw))/(1-betaa*xiw*(pitildew/piw)^(1+lambdaw/(1-lambdaw)));
        %whalohalo=((1-xiw*(pitildew/piw)^(1/(1-lambdaw)))/(1-xiw))^(1-lambdaw)...
        %    /((1-xiw*(pitildew/piw)^(lambdaw*(1+sigmaL)/(1-lambdaw)))/(1-xiw))^((1-lambdaw)/(lambdaw*(1+sigmaL)));
        smallh=whalo^(lambdaw/(1-lambdaw))*H;


        %capacity utilization parameter
        sigmab=rkbar/pinvest;
        aofu=0.5*sigmab*sigmaa+sigmab*(1-sigmaa)+sigmab*(sigmaa/2-1);

        %Creating Pitilde/Pi ratio variables
        pitildedpi=pitilded/pibar;
        pitildexpi=pitildex/pix;
        pitildemcpi=pitildemc/pimc;
        pitildemipi=pitildemi/pimi;
        pitildemxpi=pitildemx/pimx;
        pitildewpi=pitildew/piw;


        %Data series steady states
        data_pid=0; %400*log(pibar);
        data_R=400*(R-1);
        data_wdiff=0 ;
        data_cdiff=0;
        data_idiff=0;
        data_qdiff=0;
        data_Hdiff=0;
        data_ydiff=0;
        data_xdiff=0;
        data_impdiff=0;
        data_pic=0; %400*log(pic);
        data_pii=0; %400*log(pii);
        data_ystardiff=0;
        data_pistar=400*log(pistar);
        data_Rstar=400*(Rstar-1);
        data_ndiff=0;
        %data_unempdiff=0;
        %data_spreaddiff=0;
        data_spread=400*spread;
        data_gdiff=0;
        ylessd  =ybar - mu*normcdf(z-sigma,0,1)*Rk*pkprime*kbar/(pibar*muzplus);

        nxdivGDP = 0;
        Slev=.7028;
        data_Hhat=1;
        data_spreaddiff=0;
        GDP12q_diff=12*100*log(muzplus);
        GDP4q_diff=4*100*log(muzplus);
        pic4q=data_pic;
        
        %data_sdev=0;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %constructing the vector ys for Dynare%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        ys =[ log(mc); log(rkbar); log(wbar);  log(Rf); log(px); log(pix); log(k);
            log(Rnustar); log(H); log(pibar); log(mcx); log(pimc); log(mcmc);
            log(pimx); log(mcmx); log(pimi); log(mcmi); log(pic); log(pmc); log(pii); log(pmi);
            log(kbar);  log(i);  log(psizplus);  log(c);
            log(pkprime);  log(u);  log(pinvest);  %log(s); 
            a;  log(q); log(x); log(pc); log(R);
            log(ybar); log(Rx); log(pmx); log(Rk); Stilde; Sprimetilde;
            log(pitildedpi); log(Kd); log(Fd); log(phalo);
            log(pitildexpi); log(Kx); log(Fx); log(phalox);
            log(pitildemcpi); log(Kmc); log(Fmc); log(phalomc); log(cm);
            log(pitildemipi); log(Kmi); log(Fmi); log(phalomi); log(im);
            log(pitildemxpi); log(Kmx); log(Fmx); log(phalomx); log(xm);
            
            log(pitildewpi); log(Kw); log(Fw); log(whalo); %log(whalohalo); 
            log(smallh);
            
            aofu; log(piw);
          
            log(omegabar); log(n); log(spread); log(Fomegabar);
            log(muzplus); log(epsilon);  log(mupsi);  log(pitarget);
            log(Upsilon);  log(zetac);  log(Rstar); phitilde;  log(tauy); log(zetah);
            log(ystar); log(pistar);  log(g);  log(muz);
            log(taud); log(taux); log(taumc); log(taumi); log(taumx);
            log(sigma); log(gammaa);
            data_pid; data_R; data_wdiff; data_cdiff; data_idiff;
            data_qdiff; data_Hdiff; data_ydiff; data_xdiff; data_impdiff;
            data_pic; data_pii; data_ystardiff; data_pistar; data_Rstar;
            data_ndiff;%data_unempdiff; 
            data_spread; data_gdiff; log(ylessd); nxdivGDP; log(Slev);
            log(data_Hhat); data_spreaddiff; GDP12q_diff;
            GDP4q_diff; pic4q; log(pkprime*kbar-n); a/ylessd];
   

        %Dynare Switch indicating good steady state
        check=0;

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %Checking the steady states%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%

        if check==0,
            vv=[rkbar,wbar,Rnustar,Rf,Rx,R,mc,mcx,mcmc,mcmi,mcmx,pibar,pix,pic,pii, ...
                pimc,pimi,pc,px,pinvest,pmx,pmc,pmi,pkprime,kbar,sigmab,H,q,i,c,x, ...
                psizplus,ybar,ystar,cm,xm,im,n,omegabar,spread,delta,AL,gammaa,varphi];

            if min(vv) < 0
                %disp('fatal (steadystate) steady state not admissible, negative values where they are not allowed');
                ys=ys*1e+2; %passing a stupid steady state to dynare
                check=1; %indicating this is a stupid steady state; Results in a low likelihood.
            end

        end
    end

    catch
        ys=ones(size(M_.endo_names,1))*1e+2;
        check=1;
        %disp('Panic: Fsolve has a problem to terminate or other steady state calcs fail.');
        fsolve_tries=fsolve_tries+1;
   end

end


if fsolve_tries-1==max_fsolve_tries && check==1,
    disp('I give up now and return no steady state.');
end


if fsolve_tries>2 && check==0,
    %disp('Yeah, I found a steady state based on random initial guess, dude');
end



% saving s.s. values that are used as params in static file
num_para=size(M_.params,1);
for idxxx=1:1:num_para
    define_string=['M_.params(idxxx)=',M_.param_names(idxxx,:),';'];
    evalc(define_string);
end

do_steady_print=0; 
if do_steady_print==1 && check==0,
    fprintf(' \n\nsteady state properties of open economy model \n\n')
    fprintf(' Rf = %6.4f, Rx = %6.4f, R = %6.4f, mc = %6.4f, mcx = %6.4f, mcmc = %6.4f, pibar = %6.4f\n',Rf,Rx,R,1/lambdad,mcx,mcmc,pibar)
    fprintf(' pix = %6.4f, pic = %6.4f, pii = %6.4f, pimc = %6.4f, pimi = %6.4f, pc = %6.4f, px = %6.4f\n',pix,pic,pii,pimc,pimi,pc,px)
    fprintf(' pinvest = %6.4f, pmx = %6.4f, pmc = %6.4f,pkprime = %6.4f, kbar = %6.4f, sigmab = %6.4f\n',pinvest,pmx,pmc,pkprime,kbar,sigmab)
    fprintf(' q*pc = %6.4f, H = %6.4f, q = %6.4f, i = %6.4f, c = %6.4f, x = %6.4f \n',q*pc,H,q,i,c,x)
    fprintf(' a = %6.4f, psizplus = %6.4f, ybar = %6.4f, ystar = %6.4f, rkbar = %6.4f, Rk = %6.4f\n',a,psizplus,ybar,ystar,rkbar,Rk)
    fprintf(' ky = %6.4f, cy = %6.4f, iy = %6.4f, xy = %6.4f, yystar = %6.4f, real R = %6.4f\n\n',kbar/ybar,c/ybar,i/ybar,x/ybar,ybar/ystar,R/pibar)
    
    fprintf(' ****************************\n');
    fprintf(' k*Pk/y = %6.4f, c*Pc/y = %6.4f, i*Pinvest/y = %6.4f, x*Px/y = %6.4f ,n/(pkprime*k) = %6.4f\n\n',kbar*pkprime/ybar,c*pc/ybar,i*pinvest/ybar,x*px/ybar,n/(pkprime*kbar))
    fprintf(' x*capPx/(y*P) = %6.4f',px*q*pc*x/(ybar*s));
    
    fprintf(' \n model parameter values\n\n')
    fprintf(' mupsi = %6.4f, alphaa = %6.4f, pibar = %6.4f, nustar = %6.4f, nuf = %6.4f, nux = %6.4f\n',mupsi,alphaa,pibar,nustar,nuf,nux)
    fprintf(' lambdamx = %6.4f, lambdax = %6.4f, lambdad = %6.4f, omegax = %6.4f, etax = %6.4f, etag = %6.4f\n',lambdamx,lambdax,lambdad,omegax,etax,etag)
    fprintf(' AL = %6.4f, sigmaL = %6.4f, Upsilon = %6.4f, epsilon = %6.4f\n',AL,sigmaL,Upsilon,epsilon)
    fprintf(' zetac = %6.4f, zetah = %6.4f, sigmaa = %6.4f\n',zetac,zetah,sigmaa)
    fprintf(' kappaw = %6.4f, varkappaw = %6.4f, thetaw = %6.4f\n\n',kappaw,varkappaw,thetaw)

    fprintf('\n*** Financial frictions block of s.s. values below *** \n')
    fprintf('sigma = %6.4f, omegabar = %6.4f, n = %6.4f, spread = %6.4f \n\n', sigma,omegabar,n,spread );

    fprintf('*** Financial frictions block parameter values below *** \n')
    fprintf('sigmaa = %6.4f, Fomegabar = %6.4f, mu = %6.4f, we = %6.4f, gammaa = %6.4f\n\n', sigmaa,Fomegabar,mu,we,gammaa)


    varnames=[  ' ydiffU            '
        ' pidU              '
        ' RU                '
        ' wdiffU            '
        ' cdiffU            '
        ' idiffU            '
        ' qdiffU            '
        ' HdiffU            '
        ' xdiffU            '
        ' impdiffU          '
        ' picU              '
        ' piiU              '
        ' ystardiffU        '
        ' pistarU           '
        ' RstarU            '];

    means=[  100*(muzplus-1)
        400*log(pibar)
        400*log(R)
        100*(muzplus-1)
        100*(muzplus-1)
        100*(muzplus*mupsi-1)
        0
        0
        100*(muzplus-1)
        100*(muzplus-1)
        400*log(pic)
        400*log(pii)
        100*(muzplus-1)
        400*log(pistar)
        400*log(Rstar)];

    %disp('   ');
    %disp('Model means');
    %disp([varnames, num2str(means)]);
    %pause
    keyboard
end
 
