MLE program could not find feasible values

[Jasmine Xu]

Hi,

I am writing a MLE program to estimate parameters in a model following Camerer, Ho and Chong 2002 p151.
Here is my code:

program define sblmixture
    args todo b logl
    tempvar f1 f2 ff1 ff2 pp
    
    tempname alp alp_hat lamda_s lamda_bl phi_bl
    
    mleval `alp' = `b', eq(1) scalar
    mleval `alp_hat' = `b', eq(2)  scalar
    mleval `lamda_s' = `b', eq(3)  scalar
    mleval `lamda_bl'=`b', eq(4) scalar
    mleval `phi_bl'=`b', eq(5) scalar

quietly{
    gen double `f1'=.
    gen double `f2'=.
    gen double `ff1'=.
    gen double `ff2'=.
    gen double `pp'=.
    
    forval i = 1/3{
        replace s_p`i' = .
        replace s_a`i' = .
        replace s_s`i' = .
        replace bl_p`i' = .
        replace bl_s`i' = .
        replace bl_a`i' = .
    }
    
    
    *BL
    replace N=.
    by pcode: replace N = 1 if _n == 1
    by pcode: replace N = `phi_bl' * N[_n-1] + 1 if N ==.
    
    forval i = 1/3{
        by pcode: replace bl_s`i' = sigfrac`i' if _n==1
        by pcode: replace bl_s`i' = (`phi_bl'* N[_n-1]*bl_s`i'[_n-1]+sigfrac`i')/N if bl_s`i' == .
    }
    replace bl_a1=20*bl_s1+21*bl_s2+6*bl_s3
    replace bl_a2=12*bl_s1+24*bl_s2+18*bl_s3
    replace bl_a3=0*bl_s1+14*bl_s2+30*bl_s3
    
    replace bl_p1 = exp(`lamda_bl'*(bl_a1-bl_a2))/(exp(`lamda_bl'*(bl_a1-bl_a2))+1+exp(`lamda_bl'*(bl_a3-bl_a2)))
    replace bl_p2 = 1/(exp(`lamda_bl'*(bl_a1-bl_a2))+1+exp(`lamda_bl'*(bl_a3-bl_a2)))
    replace bl_p3 = exp(`lamda_bl'*(bl_a3-bl_a2))/(exp(`lamda_bl'*(bl_a1-bl_a2))+1+exp(`lamda_bl'*(bl_a3-bl_a2)))
    
    
    * S
    replace s_s1 = `alp_hat'*s_p1+(1-`alp_hat')*bl_p1
    replace s_s2 = `alp_hat'*s_p2+(1-`alp_hat')*bl_p2
    replace s_s3 = `alp_hat'*s_p3+(1-`alp_hat')*bl_p3
    
    replace s_a1 = 20*s_s1+21*s_s2+6*s_s3
    replace s_a1 = 12*s_s1+24*s_s2+18*s_s3
    replace s_a1 = 0*s_s1+14*s_s2+30*s_s3
    
    replace s_p1 = exp(`lamda_s'*(s_a1-s_a2))/(exp(`lamda_s'*(s_a1-s_a2))+1+exp(`lamda_s'*(s_a3-s_a2)))
    replace s_p2 = 1/(exp(`lamda_s'*(s_a1-s_a2))+1+exp(`lamda_s'*(s_a3-s_a2)))
    replace s_p3 = exp(`lamda_s'*(s_a3-s_a2))/(exp(`lamda_s'*(s_a1-s_a2))+1+exp(`lamda_s'*(s_a3-s_a2)))
    
    
    * Likelihood
    replace `f1' = s_p1*y1 + s_p2*y2 + s_p3*y3
    replace `f2' = bl_p1*y1 + bl_p2*y2 + bl_p3*y3
    
    by pcode: replace `ff1'=exp(sum(ln(max(`f1',1e-12))))
    by pcode: replace `ff2'=exp(sum(ln(max(`f2',1e-12))))
    
    
    replace `pp'=`alp'*`ff1'+(1-`alp')*`ff2'
    replace `pp'=. if lastobs~=1
    
    mlsum `logl'=ln(`pp') if lastobs==1
}    
end

mat start = (0.4, 0.4,  0.6, 0.6, 0.4)
ml model d0 sblmixture /alp /alp_hat /lamda_s /lamda_bl /phi_bl,tech(dfp 2 bfgs 5 nr 3)  
ml init start, copy
ml max, diff search(norescale)

When use -ml check- afrer -ml max-, I have error "could not find feasible values''. I suspect one reason could be that the line starting -replace s_a1- (also s_a2 and s_a3) because it has the probability s_p1 in it, and s_p1 is calculated later using s_a1.
But I don't know how to revise the program to accommodation this feature of the model.

Thanks a lot for your help.

[Jasmine Xu]

Another paper says because of the way that s_s1 and s_p1 are constructed "it becomes necessary of calculate a fixed point as part of fitting the model." I wonder how I can do it in Stata. Is there any other algorithm available for this type of estimation?