function [k_ideal,C] = plug_MSE_twic(lambda,b0h,sigma,maxk,stepk)

MSE(1) = inf;

for kk=0:maxk
    
    k = kk*stepk;           
    VAR = (sigma^2)*sum((1-(1-lambda).^(k+1)).^2);    
    BIAS2= sum((b0h.^2).*(1-lambda).^(2*k+2));   
    MSE(kk+2) = VAR + BIAS2;
    
    if MSE(kk+2) > MSE(kk+1)
        break;
    end
        
end
    
[C,k_twic] = min(MSE(:));

k_ideal = stepk*(k_twic-2);