Problem in transition matrix

[khoilx]

Hello Prof. Maih, I'm facing the problem of retcode "3" (Problem in transition matrix) because I declared an alien list to calculate the transitional probability (tp_1_2) from an outside function (that calculate the mvncdf of a vector) (e.g. ! s1_tp_1_2 =tp_func(...)). Could you please help me to deal with this problem in RISE?

[jmaih]

Hi Khoi,

Are you able to evaluate the function in matlab? Do you pass the alien function correctly? e.g. m=set(m, 'alien_list',{'yourFunctionName'});

J.

On Fri, Oct 11, 2024 at 1:39 AM khoilx @.***> wrote:

Hello Prof. Maih, I'm facing the problem of retcode "3" (Problem in transition matrix) because I declared an alien list to calculate the transitional probability (tp_1_2) from an outside function (that calculate the mvncdf of a vector) (e.g. ! s1_tp_1_2 =tp_func(...)). Could you please help me to deal with this problem in RISE?

— Reply to this email directly, view it on GitHub https://github.com/jmaih/RISE_toolbox/issues/189, or unsubscribe https://github.com/notifications/unsubscribe-auth/AATKBT3XWSHIUBCXSDIFB6DZ24FZXAVCNFSM6AAAAABPX2BP4GVHI2DSMVQWIX3LMV43ASLTON2WKOZSGU4DAMBWGI3DMNQ . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[khoilx]

Hello Prof. Maih, Yes, I pass the alien function correctly (I declared as: m=rise('cmr_mod','steady_state_file','model_ss',... 'rise_flags',model_settings,... 'alien_list','tp_func2');). I'm sending you the 'tp_func2.m' function file. Could you please check it for me?

function y=tp_func2(z1,z2,z3,a1,a2,a3,xmin,xmax,ymin,ymax,zmin,zmax,wmin,wmax,varargin)

if ~isreal(z1) && ~isreal(z2) && ~isreal(z3)

iz1=imag(z1);

iz2 = imag(z2);

iz3 = imag(z3);

z1=real(z1);

z2=real(z2);

z3=real(z3);

if abs(iz1)>1e-12 && abs(iz2)>1e-12 && abs(iz3)>1e-12 

    y=nan;

    return

end

end

Pev_sq = [a1 a2; a2 a3];

Pvv = [1, z3; z3, 1];

A = [z1 0; 0 z2];

try % Try solving Sigw = ASigwA'+Pvv for Sigw. Sigw = dlyap(A,Pvv); catch % If A is not stable, throw error and skip all calculations. warning('A is not stable. Assigning 0 to likelihood'); Sigw = eye(2)*.1e10; %replace s.Sigw by garbage %lik = 0; %replace s.lik by 0 end

%compute the covariance Pvve Pvve = Pvv - Pev_sq;

%compute covariance matrix for 4d multivariate normal in the calculation of time-varying transition Sigma_mvn = [Sigw, SigwA'; ASigw, Pvve + ASigwA'];

if isempty(varargin)

y=main();

else

if length(varargin)==2

    if ~strcmp(varargin{1},'diff')

        error('Wrong call to function')

    end

    switch varargin{2}

        case 1

            y=linhbin_z();

        case {2,3,4,5}

            error('Higher-order derivatives not implemented')

    end

else

    error('Wrong call to function')

end

end

function y=main()

    y=mvncdf([xmin,ymin,zmin,wmin],[xmax,ymax,zmax,wmax],[0 0 0 0],Sigma_mvn)/mvncdf([zmin,wmin],[zmax,wmax],[0 0],Sigw);

end

function fz=linhbin_z()

    fz=(mvnpdf([xmax,ymax,zmax,wmax],[0 0 0 0],Sigma_mvn)-mvnpdf([xmin,ymin,zmin,wmin],[0 0 0 0],Sigma_mvn))/...
        (mvnpdf([zmax,wmax],[0 0],Sigw)-mvnpdf([zmin,wmin],[0 0],Sigw));

end

end

And in .rs file, I define the transition probabilities as, e.g. ! s1_tp_1_2 =1-tp_func2(alpha1,alpha2,rho_v,pev11_0,pev12_1,pev22_0,-inf,b1_0,-inf,b2_0,-inf,tau1,-inf,tau2)

Thank you for your help, Khoi

[jmaih]

Hi Khoi,

I guess you are able to evaluate your function directly in Matlab. If so, one thing you may to look at is the implication of the steady state of the model for the evaluation of your function. You can always place a keyboard inside your function so that when RISE calls it, you will be able to easily see what inputs the function receives and whether it returns a finite value.

Please try this and let me know,

J-

On Fri, Oct 11, 2024 at 4:52 PM khoilx @.***> wrote:

Hello Prof. Maih, Yes, I pass the alien function correctly (I declared as: m=rise('cmr_mod','steady_state_file','model_ss',... 'rise_flags',model_settings,... 'alien_list','tp_func2');). I'm sending you the 'tp_func2.m' function file. Could you please check it for me?

function y=tp_func2(z1,z2,z3,a1,a2,a3,xmin,xmax,ymin,ymax,zmin,zmax,wmin,wmax,varargin)

if ~isreal(z1) && ~isreal(z2) && ~isreal(z3)

iz1=imag(z1);

iz2 = imag(z2);

iz3 = imag(z3);

z1=real(z1);

z2=real(z2);

z3=real(z3);

if abs(iz1)>1e-12 && abs(iz2)>1e-12 && abs(iz3)>1e-12

y=nan;

return

end

end

Pev_sq = [a1 a2; a2 a3];

Pvv = [1, z3; z3, 1];

A = [z1 0; 0 z2];

try % Try solving Sigw = ASigwA'+Pvv for Sigw. Sigw = dlyap(A,Pvv); catch % If A is not stable, throw error and skip all calculations. warning('A is not stable. Assigning 0 to likelihood'); Sigw = eye(2)*.1e10; %replace s.Sigw by garbage %lik = 0; %replace s.lik by 0 end

%compute the covariance Pvve Pvve = Pvv - Pev_sq;

%compute covariance matrix for 4d multivariate normal in the calculation of time-varying transition Sigma_mvn = [Sigw, SigwA'; ASigw, Pvve + ASigwA'];

if isempty(varargin)

y=main();

else

if length(varargin)==2

if ~strcmp(varargin{1},'diff')

    error('Wrong call to function')

end

switch varargin{2}

    case 1

        y=linhbin_z();

    case {2,3,4,5}

        error('Higher-order derivatives not implemented')

end

else

error('Wrong call to function')

end

end

function y=main()

y=mvncdf([xmin,ymin,zmin,wmin],[xmax,ymax,zmax,wmax],[0 0 0 0],Sigma_mvn)/mvncdf([zmin,wmin],[zmax,wmax],[0 0],Sigw);

end

function fz=linhbin_z()

fz=(mvnpdf([xmax,ymax,zmax,wmax],[0 0 0 0],Sigma_mvn)-mvnpdf([xmin,ymin,zmin,wmin],[0 0 0 0],Sigma_mvn))/...
    (mvnpdf([zmax,wmax],[0 0],Sigw)-mvnpdf([zmin,wmin],[0 0],Sigw));

end

end

And in .rs file, I define the transition probabilities as, e.g. ! s1_tp_1_2 =1-tp_func2(alpha1,alpha2,rho_v,pev11_0,pev12_1,pev22_0,-inf,b1_0,-inf,b2_0,-inf,tau1,-inf,tau2)

Thank you for your help, Khoi

— Reply to this email directly, view it on GitHub https://github.com/jmaih/RISE_toolbox/issues/189#issuecomment-2407581864, or unsubscribe https://github.com/notifications/unsubscribe-auth/AATKBT6BGFXVEZOB4GKN7RDZ27Q2LAVCNFSM6AAAAABPX2BP4GVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMBXGU4DCOBWGQ . You are receiving this because you commented.Message ID: @.***>

[khoilx]

Hello Prof. Maih, I realize that I can only solve the model if I use the .m function with only scalar (the normpdf of a value) to calculate the transitional probabilities, it doesn't work (as retcode shows "3" (Problem in Transition Matrix)) when I declare command of multiple variables cdf (mvdcdf). Do you have any suggestion?

[jmaih]

Hi Khoi,

Whatever you do inside the alien function is irrelevant to RISE as long as it returns a scalar and the scalar is finite. Is this the case?

On Sat, Oct 12, 2024 at 5:21 AM khoilx @.***> wrote:

Hello Prof. Maih, I realize that I can only solve the model if I use the .m function with only scalar (the normpdf of a value) to calculate the transitional probabilities, it doesn't work (as recode shows "3" (Problem in Transition Matrix)) when I declare command of multiple variables cdf (mvdcdf). Do you have any suggestion?

— Reply to this email directly, view it on GitHub https://github.com/jmaih/RISE_toolbox/issues/189#issuecomment-2408329495, or unsubscribe https://github.com/notifications/unsubscribe-auth/AATKBT2WBDEJM22HULFLA4TZ3CIU7AVCNFSM6AAAAABPX2BP4GVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMBYGMZDSNBZGU . You are receiving this because you commented.Message ID: @.***>

[khoilx]

So what is the way to use the "mvncdf" to calculate the cdf of a vector in the alien function? The output of "mvncdf" is scalar, but its input needs mean vector and covariance matrix.

[jmaih]

RISE is not going to give you scalar inputs and expects that you return a scalar output. You are free to list the covariance elements as inputs but RISE doesn't need to know how those inputs are going to be used.

On Sat, Oct 12, 2024, 10:52 khoilx @.***> wrote:

So what is the way to use the "mvncdf" to calculate the cdf of a vector in the alien function? The output of "mvncdf" is scalar, but its input needs mean vector and covariance matrix.

— Reply to this email directly, view it on GitHub https://github.com/jmaih/RISE_toolbox/issues/189#issuecomment-2408590957, or unsubscribe https://github.com/notifications/unsubscribe-auth/AATKBT6LXTE6MCAUHDFRHRDZ3EZSFAVCNFSM6AAAAABPX2BP4GVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMBYGU4TAOJVG4 . You are receiving this because you commented.Message ID: @.***>

[jmaih]

RISE will give you scalar inputs, I meant to write.

On Sat, Oct 12, 2024, 11:02 Junior Maih @.***> wrote:

RISE is not going to give you scalar inputs and expects that you return a scalar output. You are free to list the covariance elements as inputs but RISE doesn't need to know how those inputs are going to be used.

On Sat, Oct 12, 2024, 10:52 khoilx @.***> wrote:

So what is the way to use the "mvncdf" to calculate the cdf of a vector in the alien function? The output of "mvncdf" is scalar, but its input needs mean vector and covariance matrix.

— Reply to this email directly, view it on GitHub https://github.com/jmaih/RISE_toolbox/issues/189#issuecomment-2408590957, or unsubscribe https://github.com/notifications/unsubscribe-auth/AATKBT6LXTE6MCAUHDFRHRDZ3EZSFAVCNFSM6AAAAABPX2BP4GVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMBYGU4TAOJVG4 . You are receiving this because you commented.Message ID: @.***>

[khoilx]

Dear Prof. Maih, Yes, I used scalar inputs and arrange it into matrix in the alien function to calculate the multivariate normal distribution value. Could you please check if my alien function is declared correctly?

`function y=tp_func2(z1,z2,z3,a1,a2,a3,xmin,xmax,ymin,ymax,zmin,zmax,wmin,wmax,varargin)

if ~isreal(z1) && ~isreal(z2) && ~isreal(z3)

iz1=imag(z1);

iz2 = imag(z2);

iz3 = imag(z3);

z1=real(z1);

z2=real(z2);

z3=real(z3);

if abs(iz1)>1e-12 && abs(iz2)>1e-12 && abs(iz3)>1e-12

y=nan;

return

end end

Pev_sq = [a1 a2; a2 a3];

Pvv = [1, z3; z3, 1];

A = [z1 0; 0 z2];

try % Try solving Sigw = ASigwA'+Pvv for Sigw. Sigw = dlyap(A,Pvv); catch % If A is not stable, throw error and skip all calculations. warning('A is not stable. Assigning 0 to likelihood'); Sigw = eye(2)*.1e10; %replace s.Sigw by garbage %lik = 0; %replace s.lik by 0 end

%compute the covariance Pvve Pvve = Pvv - Pev_sq;

%compute covariance matrix for 4d multivariate normal in the calculation of time-varying transition Sigma_mvn = [Sigw, SigwA'; ASigw, Pvve + ASigwA'];

if isempty(varargin)

y=main(); else

if length(varargin)==2

if ~strcmp(varargin{1},'diff')

    error('Wrong call to function')

end

switch varargin{2}

    case 1

        y=linhbin_z();

    case {2,3,4,5}

        error('Higher-order derivatives not implemented')

end

else

error('Wrong call to function')

end end

function y=main()

y=mvncdf([xmin,ymin,zmin,wmin],[xmax,ymax,zmax,wmax],[0 0 0 0],Sigma_mvn)/mvncdf([zmin,wmin],[zmax,wmax],[0 0],Sigw);

end

function fz=linhbin_z()

fz=(mvnpdf([xmax,ymax,zmax,wmax],[0 0 0 0],Sigma_mvn)-mvnpdf([xmin,ymin,zmin,wmin],[0 0 0 0],Sigma_mvn))/...
    (mvnpdf([zmax,wmax],[0 0],Sigw)-mvnpdf([zmin,wmin],[0 0],Sigw));

end` And the transitional probability is defined in .rs file as following: ! s1_tp_1_2 =1-tp_func2(alpha1,alpha2,rho_v,pev11_0,pev12_1,pev22_0,-inf,b1_0,-inf,b2_0,-inf,tau1,-inf,tau2)

[khoilx]

The function uses the scalar input (e.g. 1-tp_func2(0.5, 0.5, 0.5, 0.5, 0.5, 0.5, -inf, 0, -inf, 0, -inf, 0, -inf, 0)), and give a scalar result (e.g. 0.25).

[jmaih]

1. if 0 and -inf are not going to change, there is no point in writing them as inputs. You can just define them locally inside the function
2. It is not about whether the function runs correctly when you call it with the inputs 0.5, 0.5, ....
3. What matters is when RISE calls the function, particularly the inputs RISE calls the function with. It could well be the case that RISE calls the function with extreme values such that the returned output is not finite.
4. This is why I suggested that you run RISE putting a break inside the function so that you can inspect the inputs received and check whether the output computed is finite and is between 0 and 1

On Sat, Oct 12, 2024 at 5:26 PM khoilx @.***> wrote:

The function uses the scalar input (e.g. 1-tp_func2(0.5, 0.5, 0.5, 0.5, 0.5, 0.5, -inf, 0, -inf, 0, -inf, 0, -inf, 0)), and give a scalar result (e.g. 0.25).

— Reply to this email directly, view it on GitHub https://github.com/jmaih/RISE_toolbox/issues/189#issuecomment-2408601031, or unsubscribe https://github.com/notifications/unsubscribe-auth/AATKBT6UQLQPY3TIZRHAZHLZ3E5SHAVCNFSM6AAAAABPX2BP4GVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDIMBYGYYDCMBTGE . You are receiving this because you commented.Message ID: @.***>