forwards = c(406.5, 416.64, 423.48, 433.84)
Ms = c(1, 4, 6, 9)

forward_px_w_carrying_costs = function(S, r, M, c=0.2/12){
    ##
    ## Eva1uate the price of a forward when there are carrying costs.
    ##
    one_plus_r = 1 + r/12
    discount = (1/one_plus_r)^M
    res = S/discount
    for( k in 0:(M-1) ){
        discount = (1/one_plus_r)^(M-k)
        res = res + c/discount
    }
    res
}

model_rhs = function(x){
    ##
    ## Eva1uate the difference between the quoted forward prices and the model price
    ##
    S = x[1]
    r = x[2]

    ## Compute the right-hand-side:
    ##
    rhs = rep(NA, length(forwards))
    for( mi in 1:length(Ms) ){
	rhs[mi] = forward_px_w_carrying_costs(S, r, Ms[mi])
    }

    rhs
}

meanSquaredError = function(x){
    ##
    ## Comute the MSE
    ##
    rhs = model_rhs(x)

    ## Compute the difference:
    ##
    err = rhs - forwards

    mean(err^2) ## return the MSE
}

## Solve our optimization problem:
##
S0 = mean(forwards)
r0 = 0.05

x0 = c(S0, r0)
res = optim(x0, meanSquaredError, control=list(abstol=1.e-3))
print(sprintf('S_opt= %f; r_opt= %f', res$par[1], res$par[2]))

rhs = model_rhs(res$par)
print('Predicted forward prices:')
print(rhs)