function run_MCMC_BSPF(ME_scale, nM, nSimRWMH, mhrun, nBurnRWMH,
           silent, nOutput, c, nNamePriorSpec, mData, nStartDate, nEndDate,
           nFileStr, nFileStrMH, vParameters, mSigmaProp, rnd_off, nSolutionRCLinear,
           fixed_grid,vSmolyakGridSpec,options_int_use,options_float_use)

    time_init_proc=time_ns()
    nPath = pwd()
    savedir = "$(pwd())/PosteriorDraws/" * nFileStrMH * "/"

    #--------------------------------------------------------------------------
    # Load prior specification
    #--------------------------------------------------------------------------
    include(nPath * "/PriorSpec/" * nNamePriorSpec * ".jl")

    #--------------------------------------------------------------------------
    # Configure data set and measurement errors
    #--------------------------------------------------------------------------

    vDate = mData[nStartDate:nEndDate,2]
    mY    = mData[nStartDate:nEndDate,[3,8,4,5]]
    # ygrowth (percent), log (gy), infl (annual %), R (annual %)
    println(" Start Date: $(vDate[1])")
    println(" End Date: $(vDate[end])")

    mSigmaMEtilde = [var(mY[:,1]) 0 0 0;
                     0 var(mY[:,2]) 0 0;
                     0 0 var(mY[:,3]) 0;
                     0 0 0 var(mY[:,4]) ];

    mSigmaME     = ME_scale*mSigmaMEtilde;

    #--------------------------------------------------------------------------
    # Configure proposal covariance matrix
    #--------------------------------------------------------------------------

    # set values in covariance matrix to zero that correspond
    # to fixed parameters
    vMask      = mPriorSpec[:,4]
    vFix       = mPriorSpec[:,5]
    vMaskInv   = ones(size(vMask)) - vMask
    mSigmaProp = mSigmaProp .* (vMaskInv*vMaskInv')

    #--------------------------------------------------------------------------
    # Initialize RWMH
    #--------------------------------------------------------------------------

    # Initial simulation of RWMH
    # initialize output matrices
    nStructParam = size(vParameters,1)
    mStructParam           = Array{Float64,2}(undef,nSimRWMH, nStructParam)
    vLogPosterior          = Array{Float64,1}(undef,nSimRWMH)
    vLogPosteriorCandidate = Array{Float64,1}(undef,nSimRWMH)
    mStructParamCandidate  = Array{Float64,2}(undef,nSimRWMH, nStructParam)

    #--------------------------------------------------------------------------
    # Prior / Model Solution / State-space Representation / Likelihood
    # at initial parameter value
    #--------------------------------------------------------------------------

    iSim = 1
    time_init_loop = time_ns()

    # evaluate prior at initial value
    Random.seed!(rnd_off+1)
    nLogPriorDens = fPriorDensity(vParameters,mPriorSpec)

    print(vParameters)

    # solve model and create state-space representation
    sParam = fCreateParamDict(vParameters)
    sParam_StateSpace = fStateSpaceForm(sParam,options_int_use,options_float_use,fixed_grid,vSmolyakGridSpec,nSolutionRCLinear)

    # evaluate likelihood function
    Random.seed!(rnd_off+iSim*10+2)
    (lik, s_up_stat, s_lag_stat, Eta_stat, vRegime, Neff) = fBSPF(mY, sParam_StateSpace, mSigmaME, ME_scale, nM);
    nLogLiki = sum(lik)

    print(nLogLiki)

    # compute objective function at initial value
    nObj          = nLogLiki + nLogPriorDens

    time_loop = signed(time_ns()-time_init_loop)/1000000000
    println("")
    println(" Post at Initial Para:  $nObj")
    println(" Elapsed time: $time_loop")


    #--------------------------------------------------------------------------
    # Initialization of Loop
    #--------------------------------------------------------------------------

    # initializations for loop
    nCount      = 0.0;
    nAccept     = 0.0;
    nAcceptRate = 0.0;
    nAlpha      = 0.0;
    mStructParam[1,:]          = vParameters'
    mStructParamCandidate[1,:] = vParameters'
    vLogPosterior[1]           = nObj;
    vLogPosteriorCandidate[1]  = nObj;

    # sqrt of proposal covariance Matrix
    mPropVarSVD   = svd(Matrix(c*mSigmaProp))
    mPropVar_U    = mPropVarSVD.U
    mPropVar_S    = mPropVarSVD.S
    mPropVar_sqrt = mPropVar_U * ( Matrix(1.0I, nStructParam,nStructParam).* (mPropVar_S.^(1/2)) )

    #--------------------------------------------------------------------------
    # RWMH Loop
    #--------------------------------------------------------------------------

    counter = 1
    if silent == 0
       println("")
       println(" Bayesian Computations: RWMH... ")
       println("")
    end
    time_init_loop=time_ns()

    for iSim = 2:nSimRWMH

        # Draw a THETA candidate
        Random.seed!(rnd_off+iSim*10+0)
        nCandidateShock = mPropVar_sqrt*Random.rand(MvNormal(zeros(nStructParam),Matrix(1.0I, nStructParam,nStructParam)))
        vParametersCandidate = mStructParam[iSim-1,:] + nCandidateShock
        # make sure parameter vector is consistent w fixed values
        vParametersCandidate = vMaskInv.*vParametersCandidate + vMask.*vFix

        nParametersCandidateValid = ( (vParametersCandidate[2]<=1)*(vParametersCandidate[2]>=0)*(vParametersCandidate[5]<=1)*(vParametersCandidate[6]<=1)*
            (vParametersCandidate[7]<=1)*(vParametersCandidate[8]<=1)*(vParametersCandidate[9]>0)*
            (vParametersCandidate[10]>0)*(vParametersCandidate[11]>0)*(vParametersCandidate[12]>0) )

        if  nParametersCandidateValid == 1

            # Evaluate prior at proposed draw
            nLogPriorDensCandidate = fPriorDensity(vParametersCandidate,mPriorSpec)

            # Solve model
            # Random.seed!(rnd_off+iSim*10+1)
            Random.seed!(rnd_off+1)
            sParam = fCreateParamDict(vParametersCandidate)
            sParam_StateSpace = fStateSpaceForm(sParam,options_int_use,options_float_use,fixed_grid,vSmolyakGridSpec,nSolutionRCLinear)

            # Evaluate likelihood function
            #println("*** Evaluate Likelihood Fcn ")
            Random.seed!(rnd_off+iSim*10+2)
            (lik, s_up_stat, s_lag_stat, Eta_stat, vRegime, Neff) = fBSPF(mY, sParam_StateSpace, mSigmaME, ME_scale, nM);
            nLogLikiCandidate = sum(lik)

            # objective function
            nObjCandidate = nLogPriorDensCandidate + nLogLikiCandidate
            nAlpha = min(1,exp(nObjCandidate-nObj))

        else # candidate of THETA outside of bounds
            nLogPriorDensCandidate = -1E20
            nLogLikiCandidate      = -1E20
            nObjCandidate          = nLogPriorDensCandidate + nLogLikiCandidate
            nAlpha                 = 0
        end

        # MH step
        vLogPosteriorCandidate[iSim] = nObjCandidate;
        Random.seed!(rnd_off+iSim*10+3)
        iRand = Random.rand(1);
        if iRand[1,1] <= nAlpha[1]
              mStructParam[iSim,:] = vParametersCandidate';
              nAccept              = nAccept + 1;
              nObj                 = nObjCandidate;
              vLogPosterior[iSim]  = nObjCandidate;
              nCount               = nCount + 1;
        else
              mStructParam[iSim,:] = mStructParam[iSim-1,:];
              vLogPosterior[iSim]  = nObj;
        end

        nAcceptRate = nAccept/iSim;
        mStructParamCandidate[iSim,:] = vParametersCandidate';

        counter = counter + 1
        if counter == nOutput
           time_loop = signed(time_ns()-time_init_loop)/1000000000
           println("")
           println(" Draw number:  $iSim")
           println(" Remaining draws:  $(nSimRWMH-iSim)")
           println(" Post at Current Para:  $nObj")
           println(" Post at Proposed Para: $nObjCandidate")
           println(" Acceptance rate: $nAcceptRate")
           println(" Elapsed time: $time_loop")
           counter = 0
           time_init_loop = time_ns()

           savename     = nFileStrMH * "_loop_LogPosteriorCandidate.csv"
           CSV.write(savedir * savename,  DataFrame(LogPost = vLogPosteriorCandidate[iSim-nOutput+1:iSim]), writeheader=false, append=true)

           savename     = nFileStrMH * "_loop_CandidateDraws.csv"
           CSV.write(savedir * savename,  DataFrame(mStructParamCandidate[iSim-nOutput+1:iSim,:]), writeheader=false, append=true)

           savename     = nFileStrMH * "_loop_Draws.csv"
           CSV.write(savedir * savename,  DataFrame(mStructParam[iSim-nOutput+1:iSim,:]), writeheader=false, append=true)

           savename     = nFileStrMH * "_loop_LogPosterior.csv"
           CSV.write(savedir * savename,  DataFrame(LogPost = vLogPosterior[iSim-nOutput+1:iSim]), writeheader=false, append=true)

        end

    end

    #--------------------------------------------------------------------------
    # End of RWMH Loop
    #--------------------------------------------------------------------------

    mStructParam          = mStructParam[nBurnRWMH+1:end,:];
    mStructParamCandidate = mStructParamCandidate[nBurnRWMH+1:end,:];
    vLogPosterior         = vLogPosterior[nBurnRWMH+1:end];

    # output

    savename     = nFileStrMH * "_Draws.csv"
    CSV.write(savedir * savename,  DataFrame(mStructParam), writeheader=false)

    savename     = nFileStrMH * "_CandidateDraws.csv"
    CSV.write(savedir * savename,  DataFrame(mStructParamCandidate), writeheader=false)

    savename     = nFileStrMH * "_LogPosterior.csv"
    CSV.write(savedir * savename,  DataFrame(LogPost = vLogPosterior), writeheader=false)

    savename     = nFileStrMH * "_LogPosteriorCandidate.csv"
    CSV.write(savedir * savename,  DataFrame(LogPost = vLogPosteriorCandidate), writeheader=false)

    time_total=signed(time_ns()-time_init_proc)/1000000000

    println("")
    println(" Total time: $time_total")

end
