function sum = pearstrunc(mom,mom_tar,obs,obsfac,vhuge,distyp)
%This script yields a function which is the sum of the absolute differences
%between the first four moments of a truncated Pearson family distribution
%and target moments. This function is to be minimized.

mom_act=mom; %These are the moments of the non-truncated distribution (over which optimization is to occur)

if distyp == 2 %Normal distribution: fill these values in without optimizing over them
    mom_act(3) = 0;
    mom_act(4) = 3;
end

if mom_act(4)<=mom_act(3)^2+1 %if non-truncated distribution fails condition that kurtosis > skewness^2+1
    sum=vhuge; %set value objective function to big number so it leads fminsearch away from it
else %if non-truncated distribution satisfies kurt-skew condition, then draw pearsrnd (and in general actually compute sum)
    
    %mom_act: moments of the non-truncated distribution we are working with
    %mom_t]runc: moments of the truncated distribution
    %mom_tar: the moments which we wish the truncated distribution to have
    %obs: the number of observations which we wish the truncated distribution to have
    %obsfac: indicates how many more observations we draw from the non-truncated
    %distribution to arrive at the desired number of observations.
    %e.g. obsfac=1.1 => we draw 10% more observations.
    
    %Generate one non-truncated Pearson distribution
    bigobs=round(obsfac*obs);
    pears=nan(bigobs,1);
    %rng('default');
    pears=pearsrnd(mom_act(1),mom_act(2),mom_act(3),mom_act(4),bigobs,1);
    
    %Truncate
    pearstr=nan(size(pears));
    j=0;
    for i=1:size(pears,1)
        if pears(i)>2 & pears(i)<=200 %select observations we wish to keep.
            j=j+1;  %indicates the number of observations we keep
            pearstr(j) = pears(i); %Store only values above 2 and not exceeding 200
        end
    end
    if j < obs %number of truncated observations is smaller than we need
        error('Number of truncated observations is too low. Increase obsfac.')
    elseif j >= obs %there are enough (and in general, too many) observations
        pearstr = pearstr(1:obs); %select first obs observations
    end
    
    mom_act; %These are the population moments of the non-truncated Pearson family distribution
    mom_trunc=[mean(pearstr) std(pearstr) skewness(pearstr) kurtosis(pearstr)]'; %Sample moments of the truncated Pearson
    if mom_trunc(4)<=mom_trunc(3)^2+1 %if somehow truncated dist violates skew-kurt condition (it cannot be for a dist that's already generated)
        sum=vhuge;
    else
        sum = abs(mom_trunc(1)-mom_tar(1))+abs(mom_trunc(2)-mom_tar(2))+abs(mom_trunc(3)-mom_tar(3))+abs(mom_trunc(4)-mom_tar(4))
        %SUM OF ABSOLUTE DEVIATIONS       
        %This is the objective function: minimize the sum of absolute
        %deviations of the moments of the truncated distribution from the
        %target moments (i.e. the moments of the empirical truncated
        %distribution)
    end
end
end
