Transform product into Symprod implementation

111 Views Asked by At

Heyhey,

I would like to transform the following equation for $x=(x_1,\dots, x_d)$ into a efficient matlab code:

$\begin{align}f(x) = p(e^x) \prod_{k=1}^d \varphi_{l_k,i_k}(x_k) \end{align} $

whereas $p(e^x) = p(e^{x_1}, \dots, e^{x_d}) = max \{ \sum_{k=1}^d e^{x_d}-K,0 \}$ for a given $K>0$.

Right now I have the following implementation in matlab:

phi = @(x,lvl,ix) max(1-abs(2^-lvl * x- ix),0);
payoff =@(x) max(1/d * sum(x)-K,0);
f=@(x) payoff(exp(x)) * symprod(@(k) phi(x(k),levels(k),indices(k)), k, 1, d);

Sadly this yields the error: Undefined function or variable 'k'.

Error in sparsegrids>@(x)payoff(exp(x))*symprod(@(k)phi(x(k),levels(k),indices(k)),k,1,d)

whereas, I am not sure why it says that I dont have $k$ defined.

Thanks for any suggestions.


A small working example but very poorly programmed:

for i = 1:sparse_gp
    [levels, indices] = integer2gridpoint(i,d);
    if (d == 2)
        f=@(x) payoff(exp(x)) * phi(x(1),levels(1),indices(1)) * phi(x(2),levels(2),indices(2)); 
    elseif (d == 3)
        f=@(x) payoff(exp(x)) * phi(x(1),levels(1),indices(1)) * phi(x(2),levels(2),indices(2)) * phi(x(3),levels(3),indices(3)); 
    elseif (d == 4)
        f=@(x) payoff(exp(x)) * phi(x(1),levels(1),indices(1)) * phi(x(2),levels(2),indices(2)) * phi(x(3),levels(3),indices(3)) * phi(x(4),levels(4),indices(4)); 
[...]
1

There are 1 best solutions below

1
On BEST ANSWER

This one should hopefully work.

% the dimension of x
d = 3;


% let assume some lvl and ix arrays and K
levels = rand(d,1);
indices = rand(d,1);
K=1;

% define the internal function
% note that there is no division by d in the equation for p(x)!
phi = @(x,lvl,ix) max(1-abs(2^-lvl * x- ix), 0);
payoff = @(exp_x) max(1/d*sum(exp_x) - K, 0); 

f1 = @(x) payoff(exp(x)) * prod(arrayfun(phi, x, levels, indices));

f2 = @(x) payoff(exp(x))*prod(max(1-abs(2.^(-levels).*x- indices), 0));

x = rand(d,1);
disp([f1(x) f2(x)]);

% evaluate computation time
N=1e5;
tmp=[];
tic;
for i=1:N
    tmp=f1(x);
end
t1 = toc;

tic;
for i=1:N
    tmp=f2(x);
end
t2 = toc;

fprintf("f1 takes %3.2f second, and f2 takes %3.2f second.\n", t1, t2);

Note also that you can make it more general by replacing d with length(x).