adraw.c

static double adraw (struct aspacestruct *aspace)
{   /*returns the next value in the stream of guesses of A
      terminates the sequence by sending -INFINITY
    */    
    double aa,bb,mstep,astep,upper;
    int i;
    
    aspace[0].ncalls+=1.0;
    if (aspace[0].ncalls>=ngridmax)
    {
        printf("Aspace dump:\n it=%d\n ist=%d\n id=%d\n ntogenerate=%d\n ngenerated=%d\n ncalls=%d\n baseM=%f\n baseA=%f\n lastAguess=%f\n M=%f\n keep=%d\n",
        aspace->curr->it,aspace->curr->ist,aspace->curr->id,
        aspace->ntogenerate,
        aspace->ngenerated,
        aspace->ncalls,
        aspace->baseM,
        aspace->baseA,
        aspace->lastAguess,
        aspace->M,
        aspace->keep);
        mexWarnMsgTxt("Emergiency exit from adraw, possibly infinite loop! Check model specifications!");
        aspace[0].lastAguess=-INFINITY; //stop
        return aspace[0].lastAguess;
    }
    if (aspace[0].ngenerated==0)
    {   //initial search of Aguess such that M(Aguess)<=mmax
        aspace[0].keep=0;//don't keep calculations signal
        if (aspace[0].M==INFINITY)
        {   //first step
            aspace[0].lastAguess=mmax;
        }
        else
        {
            if (aspace[0].M<=mmax)
            {   //found base point!
                // what is the base point?
                aspace[0].baseA=aspace[0].lastAguess; //abscissa (x)
                aspace[0].baseM=aspace[0].M; //ordinate (y)
                aspace[0].ngenerated=1.0; //mark change of stage
                aspace[0].keep=1; //keep calculations signal
                aspace[0].lastAguess=a0; //return a0
                aspace[0].k3=0;//init
                aspace[0].M1=a0-1;//init
            }
            else
            {   //keep looking: binary search towards a0
                if (aspace[0].lastAguess-a0<TOLERANCE)
                {
                    printf("\nLast exploratory call from adraw: A=%5.8f --> M=%5.8f, |A-a0|=%5.8e\n",aspace[0].lastAguess,aspace[0].M,fabs(aspace[0].lastAguess-a0));
                    error("Could not complete initial stage in adraw()..\nSeems like M(a0)>mmax! Increase mmax!"); 
                    return -1.0;
                }
                aspace[0].lastAguess=(aspace[0].lastAguess+a0)/2;
            }       
        }
    }    
    else
    {   //generation of next Aguess with successive linear approximations of the point where M(Aguess)==mmax
        //from now on aspace[0].ngenerated is the number of Aguess generated

        aa=(aspace[0].M-aspace[0].baseM)/(aspace[0].lastAguess-aspace[0].baseA);
        //linear model to forcast the upper bound of A guesses to be generated
        //TEMPORARY switch off rescaling and forcasting
        if (aspace[0].ngenerated==1)
        //if (aa>0 && aspace[0].M>aspace[0].M1 && aspace[0].ngenerated<aspace[0].ntogenerate-10) //last points separatelly
        {   //update upper ONLY if new point is reasonable
            bb=aspace[0].baseM-aa*aspace[0].baseA;
            upper=MIN(mmax,(mmax-bb)/aa);//updated upper bound estimate
        }
        else
        {   //or computer aa,bb from limit saved on the initialization stage
            upper=aspace[0].lim2p;
            aa=(mmax-aspace[0].M)/(upper-aspace[0].lastAguess);
            bb=mmax-aa*upper;
        }    
        aspace[0].M1=aspace[0].M;//remember returned M

        //a0 case: initializations
        //TODO: re-design how ZEROCONSUMPTION case works here and in the main body of the program
        if (aspace[0].ngenerated==1 && aspace[0].k3==0) //k3!=0 means that limits are already defined in ZEROCONSUMPTION case
        {   //initializevariables after first 'keep' call
            aspace[0].ntogenerate=(double)ngridm; //standard number of point for the above credit constraint area
            //prepare for symmetric-log grid
            aspace[0].lim2p=MIN(mmax,(mmax-bb)/aa); //upper A (plain)
            aspace[0].lim3p=-bb/aa; //focal point (plain)
            if (a0<0 && a0<aspace[0].lim3p) //sym-log grid with focal point
               aspace[0].k3=MAX( floor(aspace[0].ntogenerate*(aspace[0].lim3p-a0)/(aspace[0].lim2p-a0)), 2.0);
            else
            { //usual log grid, redefine lim3
                aspace[0].lim3p=a0;
                aspace[0].k3=1.0;
            }
            //define transformed lim
            aspace[0].lim1=tr(aspace[0].curr,aspace[0].lim3p-a0);
            aspace[0].lim2=tr(aspace[0].curr,aspace[0].lim2p-aspace[0].lim3p);
            aspace[0].lim3=tr(aspace[0].curr,0);

        }
        //zeroconsumption adjustment or standard case or stop
        if (aspace[0].M<=a0-1+TOLERANCE)
        {   //signal of c1<=0, resend point already prepared in aspace[0].lastAguess
            aspace[0].keep=1;
            //ASSUMING inverted budget calculation (in main body of the program) found exact focal point, reset it
            //also, don't need points left of it
            //also, aa and bb are corrupt because both .lastAguess and .M are wrong
            aa=(a0-aspace[0].baseM)/(a0-aspace[0].baseA);
            bb=aspace[0].baseM-aa*aspace[0].baseA;
            aspace[0].lim2p=MIN(mmax,(mmax-bb)/aa); //upper A (plain)
            aspace[0].lim3p=aspace[0].lastAguess-ZEROCONSUMPTION; //focal point (plain)
            aspace[0].k3=1.0;
            aspace[0].lim1=tr(aspace[0].curr,aspace[0].lim3p-a0);
            aspace[0].lim2=tr(aspace[0].curr,aspace[0].lim2p-aspace[0].lim3p);
            aspace[0].lim3=tr(aspace[0].curr,0);
            
        }
        else if (aspace[0].M<mmax && aspace[0].ngenerated<aspace[0].ntogenerate+0)  //allow few extra points
        {   //generate next Aguess
            aspace[0].keep=1;//keep calculations signal
            //step based on log-scale over A with linear rescaling using newer approximations of max A to be generated
            if ((int)aspace[0].ngenerated<(int)aspace[0].k3-1)
                astep=-trinv(aspace[0].curr, aspace[0].lim3
                                             +(aspace[0].k3-1-aspace[0].ngenerated)
                                             *(aspace[0].lim1-aspace[0].lim3)/(aspace[0].k3-1) 
                                                              ) +aspace[0].lim3p - aspace[0].lastAguess;
            else
                astep= trinv(aspace[0].curr, aspace[0].lim3
                                             +(aspace[0].ngenerated-aspace[0].k3+1)
                                             *(aspace[0].lim2-aspace[0].lim3)/(aspace[0].ntogenerate-aspace[0].k3) 
                                                              ) +aspace[0].lim3p - aspace[0].lastAguess;
if (astep<0)
{
    printf("\nadraw call %d\n",aspace[0].ncalls);
    printf("it=%d ist=%d id=%d\n",aspace[0].curr->it,aspace[0].curr->ist,aspace[0].curr->id);
    printf("generated %d points out of %d\n",aspace[0].ngenerated,aspace[0].ntogenerate);
    printf("aa=%1.3f\n",aa);
    printf("bb=%1.3f\n",bb);
    printf("astep=%1.3f\n",astep);
    printf("upper=%1.3f\n",upper);
    printf("lim2p=%1.3f\n",aspace[0].lim2p);
    if (astep*(upper-aspace[0].lastAguess)/(aspace[0].lim2p-aspace[0].lastAguess)<0) mexWarnMsgTxt("Step < 0");
    astep=MAX(astep,1e-5);
}            

//            aspace[0].lastAguess+=astep*(upper-aspace[0].lastAguess)/(aspace[0].lim2p-aspace[0].lastAguess);
// TEMPORARILY : switch off scaling
            aspace[0].lastAguess+=astep;
            aspace[0].ngenerated+=1.0;//increment the counter
        }
        else
        {
            aspace[0].lastAguess=-INFINITY; //stop
        }    
    }
    return aspace[0].lastAguess;
}

findLowestA.jl

function findLowestA(m::Model,cashn::Float64,p::Param,it::Int)
# check if next period consumption is positive given a level of next period cash on hand
# test is: assuming the worst case, i.e. tomorrow you get the lowest possibel income,
# what is the smallest a(t) such that c(t+1) > 0?
	
	# prepare interpolation of c(t+1)
	tmpx = vb(m.M[it+1])	# get next periods cash on hand with its lower bound
	tmpy = vb(m.C[it+1])	# get next periods consumption with its lower bound

	# if implied consumption at next period cash-on-hand cashn is negative:
	if linearapprox(tmpx,tmpy,cashn) < 0
		# find cashnext such that c=p.cfloor > 0
		cashn = p.cfloor + tmpx[1] - (tmpy[1] * (tmpx[2]-tmpx[1])) / (tmpy[2] - tmpy[1])
	end
	# return implied level of end-of-period asset that makes next period cash in hand = 0
	return invcashnext(cashn,income(m.yvec[1],it+1,p),p)	# find corresponding asset level
end