Plot a function in stata help

[Madu Abuchi]

Hi there,

Please can someone help me. I am trying to plot this function in stata but it keeps returning error message:

twoway function z=(24*y^2)/(8+y^3)^2

error in expression: (24*y^2)/(8+y^3)^2
r(111);

Any help is appreciated. Thanks

[Nick Cox]

twoway function in Stata expects a function to be defined in terms of x -- always, and regardless of whether you have any variable called x. It's just a generic placeholder for what will be plotted on the x axis. You can call the result almost anything you like, so z is fine.

This is documented, e.g.

help twoway function

So change y to x everywhere. You'll probably need to add a range() option too.

[Madu Abuchi]

Hi Nick,

Thank you very much. your response is very helpful

[Anustup Kundu]

Dear Nick Cox ,

I am trying to plot the individual prelec probability weighting function (1998) for Cumulative Prospect Theory . The function is defined as -
\[ w(p)=exp(-{(-ln(p)}^{\alpha}) \] where P is the objective probability lies between 0 and 1. W(p) is the weighting function. Alpha captures the non-linear probability weighting. I want to create a graph similar to the one on the below left using my data.

 * Example generated by -dataex-. For more info, type help dataex clear input long id float alphaa 51901 1.45 31601    1 51702 1.45 13602 1.05 32001   .6 30701   .9 23702  .35 30502  .75 11802   .8 12101    . 50501   .6 20502  .35 41502  .65 23506  .45 50701 1.45 11202   .9 53201   .6 43001 1.45 31101   .6 11302  1.3 end

If anything is unclear or if you need anything else, please let know. Best, Anustup

[Nick Cox]

Instead of your equation. I see only [Main processing error]. You can use calls to twoway function if you have an equation.

[Anustup Kundu]

Hi Nick,

My apologies for it, It was showing me but here is the equation.

w(p)=exp(-{(-ln(p)}^{\alpha})

I tried using twoway function y=exp(-(-ln(1)^alphaa)), range(0 1)

but I am not getting the graph that I referred to above. Ay idea where I went wrong?

[Nick Cox]

You need a generic x in there as the x axis argument. Its default range from 0 to 1 is exactly what you want.

twoway function exp(-((-ln(x))^2))

for example produces a nice monotonic curve. See also #2 above.

[Anustup Kundu]

Hi Nick,

Many thanks!

twoway function exp(-((-ln(x))^2)) this shows similar to what I needed.

When I plotted twoway function exp(-((-ln(x))^alphaa )) I got the graph for all the ID's. How I do the following -

1. Twoway functional graph for each ID overlayed in a single graph with the mean highlighted
2. Include a 45 degree line.

Best,
Anustup

[Nick Cox]

I guess each curve requires a separate function call -- or you create a dataset with enough data points to draw smooth curves using twoway mspline.

[Anustup Kundu]

Hi Nick,

Many thanks for your help. I did the following and I was wondering if there was any elegant way of doing it? Since the value of alpha increases in 0.05 increments. It was still possible to do it manually but if there are more points, is there any solution?

twoway ///
function exp(-((-ln(x))^0.05)), lcolor(gs13) || function exp(-((-ln(x))^0.10)) , lcolor(gs13) || ///
function exp(-((-ln(x))^0.15)), lcolor(gs13) || function exp(-((-ln(x))^0.20)) , lcolor(gs13) || ///
function exp(-((-ln(x))^0.25)), lcolor(gs13) || function exp(-((-ln(x))^0.30)) , lcolor(gs13) || ///
function exp(-((-ln(x))^0.35)), lcolor(gs13) || function exp(-((-ln(x))^0.40)) , lcolor(gs13) || ///
function exp(-((-ln(x))^0.45)), lcolor(gs13) || function exp(-((-ln(x))^0.50)) , lcolor(gs13) || ///
function exp(-((-ln(x))^0.55)), lcolor(gs13) || function exp(-((-ln(x))^0.60)) , lcolor(gs13) || ///
function exp(-((-ln(x))^0.65)), lcolor(gs13) || function exp(-((-ln(x))^0.70)) , lcolor(gs13) || ///
function exp(-((-ln(x))^0.75)), lcolor(gs13) || function exp(-((-ln(x))^0.80)) , lcolor(gs13) || ///
function exp(-((-ln(x))^0.85)), lcolor(gs13) || function exp(-((-ln(x))^0.90)) , lcolor(gs13) || ///
function exp(-((-ln(x))^0.95)), lcolor(gs13) || function exp(-((-ln(x))^1.00)) , lcolor(gs13) || ///
function exp(-((-ln(x))^1.05)), lcolor(gs13) || function exp(-((-ln(x))^1.10)) , lcolor(gs13) || ///
function exp(-((-ln(x))^1.15)), lcolor(gs13) || function exp(-((-ln(x))^1.20)) , lcolor(gs13) || ///
function exp(-((-ln(x))^1.25)), lcolor(gs13) || function exp(-((-ln(x))^1.30)) , lcolor(gs13) || ///
function exp(-((-ln(x))^1.35)), lcolor(gs13) || function exp(-((-ln(x))^1.40)) , lcolor(gs13) || ///
function exp(-((-ln(x))^1.45)), lcolor(gs13) || function exp(-((-ln(x))^1.50)) , lcolor(gs13) || ///
function exp(-((-ln(x))^0.79)), lcolor(gs2) || function exp(-((-ln(x))^1)), ///
lpattern(dash) lcolor(red) legend(off) xtitle(Objective probability) ytitle( "W(P)")

[Nick Cox]

This is a bit shorter but must be issued all at once

local call

forval j = 1/30 {
    local J = `j'/20
    local call `call' function exp(-((-ln(x))^`J')), lcolor(gs13) ||
}

twoway `call' || function exp(-((-ln(x))^0.79)), lcolor(gs2) || function exp(-((-ln(x))^1)), ///
lpattern(dash) lcolor(red) legend(off) xtitle(Objective probability) ytitle( "W(P)")

Note the loop over integers and then the division that maps 1(1)30 to 0.05(0.05)1.50. Loops over non-integers can sometimes give puzzling results. More at https://journals.sagepub.com/doi/pdf...867X1001000115

[Anustup Kundu]

Hi Nick,

Many thanks! This looks much cleaner and elegant. So as along as I loop over integers, I shouldn't expect puzzling results, right?

[Nick Cox]

Not in this context.

[Anustup Kundu]

Thanks Nick!

[Anustup Kundu]

Hi Nick,

Many thanks for your help earlier. I tried plotting the value function of for Cumulative Prospect Theory which is defined as --
\[ v(x)=x^\sigma \mathrm{for}\quad x>0\\
v(x)= - \lambda (- x)^\sigma \mathrm{for} \quad x<0
\] or if you can't view the equation then -

v(x) = x^alpha for x>0. Range for x is 0 to 1700
or, v(x) =(- lambda)*(- x)^alpha for x<0. range for x is -21 to 0. Where alpha is same as #4. Lambda is loss aversion parameter which has value 0 to 12 (continuous).

I managed to create it for a particular alpha and particular lambda (Code below). Now my questions are the following -

1. How to generalize for all alpha and all lambda. The alpha part I managed to do using your earlier code. Not sure how to extend to lambda (data below)
2. How to restrict the y-axis. So that lower values are also visible? (Check image 2).

Basic code

twoway function x^0.3, lpattern(solid) lcolor(blue) range (0 1700) || ///
function -2.19*((-x)^0.3), lpattern(solid) lcolor(red) range(-560 0) ///
ylabel(, nogrid) xline(0, lpattern("-###") lcolor(gs8)) ///
yline(0, lpattern("-###") lcolor(gs8)) ///
ytitle("Subjective value") xtitle("Objective outcome") ///
legend(order(1 "Gains" 2 "Losses") ring(0) pos(4))

Other code

local call
forval j = 1/30 {
local J = `j'/20
local call `call' function x^`J', lpattern(solid) lcolor(gs12*.5) range (0 1700) || function -2.5*((-x)^`J'), lpattern(solid) lcolor(gs12*.5) range(-560 0) ylabel(, nogrid) ||
}
twoway `call'|| function x^0.79, lpattern(solid) lcolor(gs1) range (0 1700) || ///
function -2.19*((-x)^0.79), lpattern(solid) lcolor(gs1) range(-560 0) ///
ylabel(, nogrid) yscale(range(-500 3000)) xline(0, lpattern("-###") lcolor(gs8)) ///
yline(0, lpattern("-###") lcolor(gs8)) ///
ytitle("Subjective value") xtitle("Objective outcome") ///
legend(order(1 "Gains" 2 "Losses") ring(0) pos(4))

* Example generated by -dataex-. For more info, type help dataex
clear
input long id float(alphaa lambda)
51901 1.45   1.98615
31601    1  7.601462
51702 1.45 .18429166
13602 1.05  2.042179
32001   .6  2.024825
30701   .9 4.2193832
23702  .35 .18429166
30502  .75 1.7020375
11802   .8 .22296667
12101    .         .
50501   .6 .11586667
20502  .35  2.066475
41502  .65  6.966796
23506  .45 2.1264708
50701 1.45 .18429166
11202   .9  .7416083
53201   .6  8.148367
43001 1.45    1.5399
31101   .6 .11586667
11302  1.3  1.540396
end