-1
$\begingroup$

I am trying to repeat an initial condition many times such that I get a desired final output. ONE final output depends on many subsequent calculations from ONE initial condition. But I want to get a specific final output. I have many given parameters

Clear["Global`*"]

tid = 50; taid = 20; tcid = 30; rho = 0.975; pe1 = 4.37; pe2 = 4.4509; qedom = 1021167; 

qeim = 5700; qcid = 529.59; el = 0.655; did = 4.22; dtot = 430932.474; iad = 445002; 

cad = 688723.8554;

qetot = qedom + qeim;

uel = 1 - el;

sigid = 0.20568; rid = 0.02958;

upid = N[Exp[Sqrt[1]*sigid]];

dwid = N[1/upid];

rf = N[Exp[rid*1]];

pid = N[(rf - dwid)/(upid - dwid)]; 

qid = N[1 - pid];

prf2 = 0.001*pe2*qetot-0.01*cad; 

npva = ((rho - rho^(taid + 1))/(1 - rho))*prf2 - 0.01*iad;  

I would encourage you NOT TO CHANGE the values of the parameters(you are free to change the names if you wish); these are from real data and heavily influence the results. The following parameter "cprd" is the INITIAL CONDITION (given the other parameters) that I want to change many times

cprd = 200; 

pin = 63.56*(cprd);(*63.56 is a FIXED constant and "pin" is the product 63.56*cprd*)

prf1[p_] = 0.001*pe1*qedom*el - 0.001*p*qcid - 0.001*pe1*qeim;

As you see, ALL THE FOLLOWING calculations depends on ONE chosen value of "pin" and if I change "cprd" all the results would change.

stsp = Table[pin*(upid^i), {i, -tid, tid}];

Dimensions[stsp] (*checking this gives a column vector*)

{101}

lensp = Length[stsp]

101

prf1[stsp];

per[i_] := prf1[stsp][[i]];

up[n1_, sig_, T_] := N[Exp[Sqrt[T/n1] sig]];

down[n1_, sig_, T_] := N[1/up[n1, sig, T]];

int[n1_, Rf_, T_] := N[Exp[Rf*(T/n1)]];

P[up_, down_, al_] := N[(al - down)/(up - down)];

Q[up_, down_, al_] := N[1 - P[up, down, al]];

mean[l_List] := Apply[Plus, l]/Length[l];

AmericanOption[p0_, n1_, sig_, T_, Rf_, exercise_Function] := 
  Module[{u = up[n1, sig, T], d = down[n1, sig, T], 
    al = int[n1, Rf, T], p, q, OpRecurse, res}, p = P[u, d, al]; 
   q = Q[u, d, al]; 
   OpRecurse[node_, level_] := 
    OpRecurse[node, level] = 
     If[level == n1, exercise[p0*d^node u^(level - node)], 
      rho*{p, q}.{OpRecurse[node, level + 1], 
          OpRecurse[node + 1, level + 1]} + 
       exercise[p0*d^node u^(level - node)]]; res = OpRecurse[0, 0]; 
   Clear[OpRecurse]; res]; 

AmericanPut[p0_, n1_, sig_, T_, Rf_] := AmericanOption[p0, n1, sig, T, 
  Rf, # &] 

AmericanPut[per[1], tcid, sigid, tcid, rid]  (*just checking*)
95832.1

tabnpvc =Table[AmericanPut[per[i], tcid, sigid, tcid, rid], {i, 1, lensp}];

e1 = 0;

newp = Table[e1, {i, lensp}, {j, lensp}];

Dimensions[newp]

{101, 101}

Table[newp[[i, Min[i + 1, lensp]]] = pid, {i, lensp}]; 


Table[newp[[i, Max[i - 1, 1]]] = qid, {i, lensp}];

    valtern = Max[#, npva] & /@ tabnpvc;

The matrix "val" is the MOST IMPORTANT. I want to REPEAT THE ITERATIONS (by changing "cprd") such that difference between the 1st and the 50th column val[[All, 1]]-val[[All, 50]] is first positive and then zero (vector contains positive numbers followed by zeros).

val = Table[e1, {i, lensp}, {j, tid + 1}];

val[[All,tid + 1]] = tabnpvc; 

Do[Table[val[[All, i]] = Max[#, npva] & /@ (prf1[stsp] + rho*newp.val[[All, i + 1]]), {i, 
    1, tid}], {tid}];

val // MatrixForm; 

Dimensions[val];

diff = val[[All, 1]] - val[[All, 50]]; (*this vector should ideally contain positive values FOLLOWED BY zeros but NOT negative*)

pos = Position[diff, _?(# < 0 &)] (*this is not desired*)
{{36}, {37}, {38}, {39}, {40}, {41}, {42}, {43}, {44}, {45}, {46},{47}}

Finally, just checking with figures

valplot[i_, j_] := {stsp[[i]], val[[All, j]][[i]]}; 

val1 = ListLinePlot[Table[valplot[i, 1], {i, 1, lensp}], 
  PlotRange -> All]; 

val2 = ListLinePlot[Table[valplot[i, 50], {i, 1, lensp}], 
  PlotRange -> All]; 

Show[val1, val2]; (*this should be such that the curves don't ever cross but converge to the value of "npva"*)

So MY DESIRED FINAL OUTPUT is that the vector "diff" should contain BOTH positive numbers and zeros. It should first be positive and then zeros (maybe "pos" should equal the null vector {}). As you see, you can keep changing the value of "cprd" until you get the desired result

I have no idea on how to do incorporate all this in a "While" or "Do" function. Please help. Thank You.

$\endgroup$
  • 2
    $\begingroup$ Can't you just encapsulate all that in a function? I'll do it for you, but it's not so hard. And then you can use While with that function. It also doesn't seem like you separated out just the parts of your code that are relevant to the question. Please don't just dump your entire notebook into the question (e.g., val // MatrixForm; is a totally useless line). The simpler your question is the more able and willing people will be to answer it. $\endgroup$ – b3m2a1 Dec 30 '16 at 6:03
  • $\begingroup$ @MB1965, is the code Replace[Do[ With[{val = OPsBigFunctionThing[crpd]}, If[FreeQ[diff@val, _?(# < 0 &)], Return[crpd]]], {crpd, 300, 1000}], Null -> $Failed] searching over 300 -1000? When I run, I get "$Failed", does that mean that the vector diff[val] is null over the range 300 to 1000? Would have to change the range then? Also, could you plz explain what you mean by "optimizing?" The code diff@OPsBigFunctionThing[-0.25] // AbsoluteTiming // First gives me 1.688220 (not your 0.741599). If I am right, it just gives the number of seconds taken to compute diff when cprd=-0.25?? Tks $\endgroup$ – Supratim Das Gupta Jan 1 '17 at 1:02
2
$\begingroup$

So the simplest way to tackle this is just to put all your code into a function. You may want to optimize it, because each iteration takes a while, but that's another story.

I did that for you:

OPsBigFunctionThing[cprd_]:=
    Module[{
        tid=50,taid=20,tcid=30,rho=0.975,
        pe1=4.37,pe2=4.4509,qedom=1021167,
        qeim=5700,qcid=529.59,el=0.655,
        did=4.22,dtot=430932.474,iad=445002,
        cad=688723.8554,qetot,
        sigid=0.20568,rid=0.02958,uel,
        upid,dwid,rf,pid,qid,prf2,npva,
        pin = 63.56*(cprd)(*63.56 is a FIXED constant and "pin" is the product 63.56*cprd*),
        prf1,stsp,lensp,
        per,up,down,int,
        P,Q,mean,AmericanOption,
        AmericanPut,tabnpvc,e1,
        newp,valtern,val
        },

        qetot=qedom+qeim;
        uel=1-el;
        upid=N[Exp[Sqrt[1]*sigid]];
        dwid=N[1/upid];
        rf=N[Exp[rid*1]];
        pid=N[(rf-dwid)/(upid-dwid)];
        qid=N[1-pid];
        prf2=0.001*pe2*qetot-0.01*cad;
        npva=((rho-rho^(taid+1))/(1-rho))*prf2-0.01*iad;

        prf1[p_] := 0.001*pe1*qedom*el - 0.001*p*qcid - 0.001*pe1*qeim;

        stsp = Table[pin*(upid^i), {i, -tid, tid}];
        lensp = Length[stsp];

        per[i_] := prf1[stsp][[i]];
        up[n1_, sig_, T_] := N[Exp[Sqrt[T/n1] sig]];
        down[n1_, sig_, T_] := N[1/up[n1, sig, T]];

        int[n1_, Rf_, T_] := N[Exp[Rf*(T/n1)]];
        P[up_, down_, al_] := N[(al - down)/(up - down)];
        Q[up_, down_, al_] := N[1 - P[up, down, al]];

        mean[l_List] := Apply[Plus, l]/Length[l];

        AmericanOption[p0_, n1_, sig_, T_, Rf_, exercise_Function] := 
            Module[{
                u = up[n1, sig, T], d = down[n1, sig, T], 
              al = int[n1, Rf, T], p, q, OpRecurse, res}, 
              p = P[u, d, al];
              q = Q[u, d, al]; 
              OpRecurse[node_, level_] := 
                OpRecurse[node, level] = 
                    If[level == n1, exercise[p0*d^node u^(level - node)], 
                        rho*{p, q}.{OpRecurse[node, level + 1], 
                        OpRecurse[node + 1, level + 1]} + 
                        exercise[p0*d^node u^(level - node)]]; 
             res = OpRecurse[0, 0]; 
        Clear[OpRecurse]; 
        res]; 

        AmericanPut[p0_, n1_, sig_, T_, Rf_] := 
            AmericanOption[p0, n1, sig, T, Rf, # &];

        tabnpvc = Table[AmericanPut[per[i], tcid, sigid, tcid, rid], {i, 1, lensp}];
        e1 = 0;

        newp = Table[e1, {i, lensp}, {j, lensp}];
        Do[newp[[i, Min[i + 1, lensp]]] = pid, {i, lensp}]; 
        Table[newp[[i, Max[i - 1, 1]]] = qid, {i, lensp}];
        valtern = Max[#, npva] & /@ tabnpvc;
        val = Table[e1, {i, lensp}, {j, tid + 1}];
        val[[All,tid + 1]] = tabnpvc; 
        Do[
            Table[
                val[[All, i]] = Max[#, npva] & /@ (prf1[stsp] + rho*newp.val[[All, i + 1]]), 
                {i, 1, tid}], 
            {tid}];
        val
]

Then we'll make your diff its own function:

diff[val_] := val[[All, 1]] - val[[All, 50]]

And finally use a Do loop to try to find a value for you:

Replace[
    Do[
        With[{val=OPsBigFunctionThing[crpd]},
            If[FreeQ[diff@val,_?(#<0&)],
                Return[crpd]
                ]
            ],
        {crpd,300,1000}],
    Null->$Failed
    ]

Note that you'll probably want to optimize:

In[32]:= diff@OPsBigFunctionThing[-0.25] // AbsoluteTiming // First

Out[32]= 0.741599

Depending what sort of range you're searching over that could be excruciating.

BTW: I've already searched 1..1000 for you so you can take over from there.

$\endgroup$
  • $\begingroup$ I don't know a better way to thank you....a Million Thanks! $\endgroup$ – Supratim Das Gupta Dec 30 '16 at 15:16
  • 1
    $\begingroup$ @SupratimDasGupta A good way of saying "thank you" is to accept the answer. You should consider browsing through your other questions and accepting the answers that were useful - you haven't done that in almost all cases. $\endgroup$ – corey979 Dec 30 '16 at 15:26
  • $\begingroup$ I would be more than happy to do it. I could not find any tab "accept" in the page. $\endgroup$ – Supratim Das Gupta Dec 31 '16 at 17:05
  • $\begingroup$ There's a little check mark to click next to the answer. $\endgroup$ – b3m2a1 Dec 31 '16 at 17:38
  • $\begingroup$ @MB1965, is the code Replace[Do[ With[{val = OPsBigFunctionThing[crpd]}, If[FreeQ[diff@val, _?(# < 0 &)], Return[crpd]]], {crpd, 300, 1000}], Null -> $Failed] searching over 300 to 1000? When I run, I get "$Failed", does that mean that the vector diff[val] is null over the range 300 to 1000? Would have to change the range then? Also, could you plz explain what you mean by "optimizing?" The code diff@OPsBigFunctionThing[-0.25] // AbsoluteTiming // First gives me 1.688220 (not your 0.741599). If I am right, it just gives the number of seconds taken to compute diff when cprd=-0.25?? Tks $\endgroup$ – Supratim Das Gupta Jan 1 '17 at 1:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.