Richardson extrapolation in BS PDE

For solving PDE numerically, especially for option which has discontinuity in
initial configuration, it is hard to use high order differention function.

In our test, we find solelly, if we use O(h^2) order approximation to the differention, we can only achieve O(h) order approximation.

If we use O(h^4) order approximation, the error is roughly O(h^2), at this time,
a Richardson extrapolation can reduce the error significantly. (As in this
situation, it is not easy to reduce the space step unboundedly)

How to export the residuals in SAS

proc reg data=hartnagl;
model ftheft = fertil labor postsec mtheft/dw;
output out=hartreg p=p r=r; /* We import this to hartreg with p is prediction
r is the OLS residual*/
run;
quit;

Robust optimizaton

If our solution are subject to error(which is
not a constant), it is meaningful to consider
the so called Robust Optimization.

Comparing to global optimization, the robust
optimization may do not find the highest point
but the most stable peak.

(MIN MAX (WORST CASE)) is the robust solution.

Note 12: Finite difference approach to solve Heston model(2)

It is important for choosing method for the time propogating.

1. Explicit schemes are easily to implement, do not suffer from oscillation problems but are only conditionally stable.

2. Implicit schemes are also oscillation-free, unconditionally stable but only first-order accurate.

3. The Crank–Nicolson scheme is (theoretically) second-order accurate but it is well known that it produce spurious oscillations and spikes near the strike price, barriers and monitoring points.

4. Hybrid method. Such as in Rannacher (1984) it uses fully implicit
Euler for the first few time steps and Crank–Nicolson after that. The scheme is stable, firstorder accurate and is oscillation-free. And use Richardson extrapolation in combination with implicit Euler (Gourlay, 1980). The resulting
scheme is stable, second-order accurate and again oscillation free.

5. Function smooth.1) Projecting the initial condition onto a set of basis functions


6. One can combine smoothing and the Rannacher method to get a stable, second-order accurate and oscillation-free finite difference scheme

Note 12: Finite difference approach to solve Heston model(1)

1) By Ito formula, it is easy to write down the PDE for
Heston model
2) The inmportant thing is the boundary condition
s->0 s-> infinity, u->0 u-> infinity
a) For European call option
we have s->0 U=0 u-> infinity U=S
u=0, we solve the PDE separatly which can be regarded a boundary
s-> infinite , Let us use U'=1(linear boundary condition)
b)For European call option
we have s->0 U=K u-> infinity U=0
u=0, max(K-S,0)
s-> infinite , Let us use U'=0(linear boundary condition)
c) we can have some other conditions
3) Numerical method for Heston model

Let us give basic introduction to the numerical procedure
1) Create Mesh
2) Define initial condition
3) update the boundary conditions for time level n+1 first at time level n
4) Propogate the state at time level n to time level n+1
5) Iterative Until the time is ended.

Note 11: How to solve parital integro differential equation(PIDE)

In BS formalism, we know V_t=LV which is a partial differential equation. However, if we included jump
as Merton model, our equation changes to V_t=LV+I(V),
I(V) means an integration term. So now it is PIDE.

1) How to deal with the integration term
* FFT evaluation
* a) Interpolate the discrete values to equal space (Or using unequal space FFT evaluation, it is costly) b) after the result, interpolate back...
The reason is that we usually use unequally spaced gird in S coordinates as near the strikes and barriers we use small mesh.
* For avoid “wrap around pollution”, use linear extrapolation to extend the Smax, Smin.

2) For each step, after getting the value of integration, we can boost the calculation using Fix point interation et al. to use the PDE in this step.

Use SAS to fit the ARIMA model with the Box-Jenkins methodology

proc arima data=test;
by id
identify var=trend scan;
run;

estimate p=2 q=1 ;
run;

forecast
lead=10 out=predicttrend printall;
run;
quit;

Note 10: Detailes for the usage of QuantLib Calendar Class

When we deal with real business time, it is no more simple time interval calculation between two different time. And one need to use the business day calculation. In principle, it is extremely simple and just needs basic level of programming. However, it is tedious for a quant working in this convention. So Let us borrow others.


Date todayDate(15, May, 1998);
Date settlementDate(17, May, 1998);
Date startDate(18,May,1998);
Date maturity(17, May, 1999);
// We use the Holidays setting for the stock exchange in NYSE
// Market { Settlement, NYSE, GovernmentBond, NERC }@USA
Calendar calendar = UnitedStates();

// Calculate the business day between
int btimeinterval=calendar.businessDaysBetween(startDate,maturity);
// Calculate the day between
DayCounter daycount = Actual365Fixed();
int timeinterval=daycount.dayCount(startDate,maturity);

std::cout << "SettlementDate = " << settlementDate << std::endl;
std::cout << "StartDate = " << startDate << std::endl;
std::cout << "Maturity = " << maturity << std::endl;
std::cout << "Total business days= "<< btimeinterval << std::endl;
std::cout << "Totaldays= "<< timeinterval << std::endl;


Date startDate01 = calendar.advance(settlementDate,Period(1,Days));

Date maturity01 = calendar.advance(settlementDate,Period(1,Years));

std::cout << "SettlementDate = " << settlementDate << std::endl;
std::cout << "StartDate = " << startDate01 << std::endl;
std::cout << "Maturity = " << maturity01 << std::endl;

Note 09: Black-Scholes Option Pricing Formula And Exotic Options

1. THE BLACK-SCHOLES OPTION PRICING FORMULA (Black-Scholes (1973))

Theorem. Consider a European call option on a stock whose current price is S. Suppose that the stock price is lognormally distributed with volatility R, that the option’s exercise price is X, that the exercise date of the option is T , and that the continuously compounded interest rate is r . Furthermore assume that the stock will pay no dividends before the option exercise date T . Then the call price is given by:

C = SN (d₁) - Xe^-rT N(d₂)‚
Where d₁=(ln(S/X)+(r+σ²/2)/T)/(σ√T), d₂=d₁-σ√T. And N() means values of the cumulative
standard normal distribution. (We can use Erf[x/1.41421 3562373095048801688724]/2+0.5 to Calculate N() in Mathematica)

Also Put-Call Parity Theorem tell us P=C+ Xe^-rT -S.

So it is ready to calculate the option if we know the above parameter.
Here is the mathematical code

snormal[x_]:=Erf[x/1.414213562373095048801688724]/2+0.5;
Clear[d1, d2, bsCall, bsPut]

d1[s_, x_, sigma_, T_, r_]:=(Log[s/x]+(r+sigma^2/2)*T)/
(sigma*Sqrt[T])

d2[s_, x_, sigma_, T_, r_]:=d1[s, x, sigma, T, r]-sigma*Sqrt[T]

bsCall[s_, x_, sigma_, T_, r_]:=
s*snormal[d1[s,x,sigma,T,r]]-x*Exp[-r*T]*snormal[d2[s,x,sigma,T,r]]

bsPut[s_,x_,sigma_,T_,r_]:=bsCall[s,x,sigma,T,r]+x*Exp[-r*T]-s

Later, I will show how to calculate this using Quantlib, the practical issue for this point is 1) How to define
the date interval 2) Inversely calculate the volatility surface 3) How to calibrate the model 4) How to work our C++ program with the Excel

(to be continue)

Test of my global evolutionary optimization program

It seemed that my program works quite well for
global optimization. It even found some new global
minimization point. If possible, I would like to
improve it further and write a paper on this subject.
(For definition of those functions, see
Some new test functions for global
optimization and performance of
repulsive particle swarm method
Sudhanshu Mishra
North-Eastern Hill University, Shillong (India))

For Hougen function, our program find the solution is about
1 1.72812 0.0868823 0.0566724 0.15568 0.868809 as
the standard deviation is 0.298914. (The best result is
Rosenbrock-Quasi-Newton method:
1 1.253031 1.190943 0.062798 0.040063 0.112453. The sum
of squares of deviations is 0.298900994).


--------------------

The Hansen function
f* = -176.541793 Global Solution=9
1 4.97648 -7.70831 -176.542
2 -7.58989 -7.70831 -176.542
3 -1.30671 -1.42513 -176.542
4 -1.30671 4.85806 -176.542
Total finding is 9
The tube holder function a
x* = -10.8723 Global Solution=2
1 1.5706 5.74051e-009 -10.8723
2 -1.5706 -8.72265e-009 -10.8723
Total finding is 2
The tube holder function b
x* = -10.8723 Global Solution=2
1 1.5706 5.75833e-009 -10.8723
2 -1.5706 9.85915e-010 -10.8723
Total finding is 2
The Holder table function
x* = -26.92 Global Solution=4
1 9.64617 9.64617 -26.9203
2 9.64617 -9.64617 -26.9203
3 -9.64617 -9.64617 -26.9203
4 -9.64617 9.64617 -26.9203
Total finding is 4
The Carrom table function
x* = -24.1568155 Global Solution=4
1 -9.64617 -9.64617 -24.1568
2 -9.64617 9.64617 -24.1568
3 9.64617 9.64617 -24.1568
4 9.64617 -9.64617 -24.1568
Total finding is 4
The Cross in tray function
x* = -2.06261218 Global Solution=4
1 1.34941 -1.34941 -2.06261
2 -1.34941 -1.34941 -2.06261
3 1.34941 1.34941 -2.06261
4 -1.34941 1.34941 -2.06261
Total finding is 4
***The Crowned cross function
x* -> 0.0 Global Solution=1
1 9.99995e-016 3.27808e-018 0.0012431
Total finding is 1
***The cross function
x* -> 0.0 Global Solution=4
1 -1.34941 -1.34941 4.84822e-009
2 1.34941 -1.34941 4.84822e-009
3 1.34941 1.34941 4.84822e-009
4 -1.34941 1.34941 4.84822e-009
Total finding is 4
*****Cross-leg table function:
x* -> -1.0 Global Solution=1
1 -6.74913e-018 9.99978e-016 -7.48396e-006
Total finding is 1
Pen holder function:
x* -> -0.96354 Global Solution=1
1 9.64617 9.64617 -0.963535
2 9.64617 -9.64617 -0.963535
3 -9.64617 9.64617 -0.963535
4 -9.64617 -9.64617 -0.963535
Total finding is 4
Bird function:
x* -> -106.764537 Global Solution=1
1 -1.58214 -3.13025 -106.765
2 4.70104 3.15294 -106.765
Total finding is 2
Modified Schaffer function #1:
x* = 0 Global Solution=1
1 3.24218e-007 -3.02404e-007 0
Total finding is 1
Modified Schaffer function #2:
x* = 0 Global Solution=1
1 -3.42333e-007 2.56129e-007 0
2 0.00139488 -0.00139456 0
3 0.00211379 -0.00211354 0
4 0.0024989 -0.00249896 0
Total finding is 133
Modified Schaffer function #3:
x* = 0.00156685(0,1.253115) Global Solution=1
1 41.0424 -41.0615 3.66026e-007
Total finding is 1
Modified Schaffer function #4:
x* = f (0, 1.253132) = 0.292579 Global Solution=1
1 -41.7439 41.7251 0.291927
Total finding is 1
Styblinski-Tang function
x* =f (-2.903534, -2.903534)=-78.332 Global Solution=1
1 -2.90353 -2.90353 -78.3323
Total finding is 1
Bukin function f6
x* =f (-10, 1)=0 Global Solution=1
1 -10.0009 1.00018 1.13323e-005
Total finding is 1
Giunta function:
x* =f (0.45834282, 0.45834282) =0.0602472184 Global Solution=1
(Wrong! f (0.45834282, 0.45834282) =0.0646388)
1 0.46732 0.46732 0.0644704
Total finding is 1
Schaffer function
x* =f (0, 0) =0 Global Solution=1
1 -2.06858e-010 -1.61139e-009 0
Total finding is 1
Press any key to continue

Note 08: An study on an example of QuantLib: EquityOption

1)Notice this example does not contain calibration procedure.
2) Result from executive file after compiling:

Option type = Put
Maturity = May 17th, 1999
Underlying price = 36
Strike = 40
Risk-free interest rate = 6.000000 %
Dividend yield = 0.000000 %
Volatility = 20.000000 %


Method European Bermudan American
Black-Scholes 3.844308 N/A N/A
Barone-Adesi/Whaley N/A N/A 4.459628
Bjerksund/Stensland N/A N/A 4.453064
Integral 3.844309 N/A N/A
Finite differences 3.844342 4.360807 4.486118
Binomial Jarrow-Rudd 3.844132 4.361174 4.486552
Binomial Cox-Ross-Rubinstein 3.843504 4.360861 4.486415
Additive equiprobabilities 3.836911 4.354455 4.480097
Binomial Trigeorgis 3.843557 4.360909 4.486461
Binomial Tian 3.844171 4.361176 4.486413
Binomial Leisen-Reimer 3.844308 4.360713 4.486076
Binomial Joshi 3.844308 4.360713 4.486076
MC (crude) 3.834522 N/A N/A
MC (Sobol) 3.844613 N/A N/A
MC (Longstaff Schwartz) N/A N/A 4.481675

Run completed in 5 s

3) With the above result, let us begin to analyze how this program works step by step:

Step 1: Include the QuantLib header and Boost timer
#define BOOST_LIB_DIAGNOSTIC #  include  #undef BOOST_LIB_DIAGNOSTIC
#include 

Step 2: Define option parameters: The important thing here
is that we need to the time as Date class. And also We need
to determine the date counter type(DayCounter class).
      Option::Type type(Option::Put);

       Date todaysDate(15, May, 1998);
       Date settlementDate(17, May, 1998);
       Settings::instance().evaluationDate() = todaysDate;
       Date maturity(17, May, 1999);
       DayCounter dayCounter = Actual365Fixed();

Step 3: Determine the exercise Dates for different type of options/
And also we need to set the parameters to the engine of the Quantlib.
The trick here is to use the Handle class to make the observable being shared.

template
class QuantLib::Handle<> And is inherited from observable.

So the definition for dividend is like:

        Handle flatDividendTS(
            boost::shared_ptr(
               new FlatForward(settlementDate, dividendYield, dayCounter)));
Step 4: define option type , eg: 
VanillaOption americanOption(stochasticProcess, payoff,
                                   americanExercise);
Step 5: Use different method to calculate the option. The only trick here is to
use the setPricingEngine to do this job with the selected. eg:
     americanOption.setPricingEngine(boost::shared_ptr(
                                new BaroneAdesiWhaleyApproximationEngine));

how can we output the SAS estimated parameters?

1)
For each SAS procedure, SAS produces a group of ODS output objects. In order to know what objects are associated with a particular proc, we use ods trace on statement right before the proc and turn the trace off right after it.

ods trace on /listing;
proc reg data=****;
model *****;
run;
quit;
ods trace off;
2) After knowing the label of object, we can
use the label to create a data set just as using
the name. Eg, one label of the Reg Output is
"Parameter Estimates", now we can try:
ods output "Parameter Estimates"=parest;
proc reg data=***;
 model ****;
run;
quit;
ods output close;


Trick:

A)
We can NOT use noprint option since ODS requires an output object.but the statement
ods listing close;
eliminates the output to appear in the output window. After the proc reg, we turn
 back the listing output back so output will appear in the output window again.

ods listing close;
ods output "Parameter Estimates"=****;
proc reg data=****;
model ****;
run;
quit;
ods output close;
ods listing;
 B)Output to an HTML file

Let's say that we want to write the output of our proc reg to an HTML file. This can be done very easily using ODS. First we specify the file name we are going to use. Then we point the ods html output to it. At the end we close the ods html output to finish writing to the HTML file. You can view procreg.html created by the following code.

filename myhtml "c:\examples\procreg.html";
ods html body=myhtml;
proc reg data=hsb25;
model write= female math;
run; quit;
ods html close;

Note 06: Dynamic Hedging in real world(I)

Introduction

1) Continuous Dilemma
a) Transition Cost( If we want to hedge at high frequency to mimic the continuous/theoretical method, we have to face it )
b) Variance of return in hand (different time, different situation)

2)Risk Management
1) Financial risk inherent to nonfinnancial business
2) Market risk incurred by the provider of finanical instrument
Micromanagement/Macromanagement for hedging
Residual Risks(Legal risks, Fraud risks, Credit risks....)

3) Gap between trader and quant
a) Explain the main conclusion in a single sentence before discuss the subject matter
b) Explain the subject matter in a single sentence
c) Reject the whole project if unable to perform the above two steps.

Note 05：Volatility Smile（II）

4）Universal Volatility Model
Interestingly，one can combine local volatility model( which regards the volatility as a
function of time and stock price) and stochastic volatility(jump diffusion) model together to construct a universal volatility model (Blacher/Lipton). More parameters, more accuracy...

5) Regime Switch Model
This model admit that the stock price must be modeled by different regimes which are different
parameters. In first sense, this model is not so realistic as even we allow smoothly (in time) swith from a regime to one regime, the market dynamics is not so simple as there is no reason to protect whether the real dynamics is governed by regimes. However, one can implement the above 1) 2) 3) 4) models as different regimes, with the inducement of more parameters, one can surely achieve higher precision.

6) Calibration by exotic option
As we have a lot of models in hand, every one seems to fit the volatility well. How can we do? Of course, directly calibrate the model using exotic option is a good idea.


Then there comes the question. Parameters V.S The Prediction


(to be continue)

How to analyze Panel data using SAS

Step 1: Sort the data properly
eg:
proc sort data=a;
by state date;
run;
Step 2: invoke the TSCSREG procedure and specify the cross section and time series variables in an ID statement. And specify the linear regression model with a MODEL statement: the dependent variable is listed first, followed by an equal sign, followed by the list of regressor variables

eg:
proc tscsreg data=a;
id state date;
model y = x1 x2;
run;

In order to aid in model specification within this class of models, the procedure provides two specification test statistics Rejection of the null hypothesis might suggest that the fixed effects model is more appropriate.
a) F statistic that tests the null hypothesis that the fixed effects parameters are all zero.
b) a Hausman m-statistic that provides information about the appropriateness of the random effects specification

Fixed effects: the models are essentially regression models with dummy variables corresponding to the specified effects

Random effects:
We include a classical error term with zero mean and a homoscedastic covariance matrix
One-Way model: we also include error term with the index of cross section
Two-Way model: We include also error term with the time change

Usually you cannot explicitly specify all the explanatory variables that affect the dependent variable. The omitted or unobservable variables are summarized in the error disturbances. The TSCSREG procedure used with the Fuller-Battese method adds the individual and time-specific random effects to the error disturbances, and the parameters are efficiently estimated using the GLS method

Step 3: The following statements are used with the TSCSREG procedure.
PROC TSCSREG options;
BY variables;
ID cross-section-id-variable time-series-id-variable;
MODEL dependent = regressor-variables / options;
label: TEST equation [,equation... ];

Note 04: My understanding to Volatility Smile and it's Models

It is well known for almost everyone in finance that volatility surface is not flat.
In my understanding, volatility smile is a result of market multi-interaction: Hedging(fear of risk) , Stochastic Volatility/Jump Diffusion , Different opinions on future trend and so on. But the main reason is that the emergence of the new jobs:Quant. The reason is simple. Before 1987, the surface is flat. After 1987, with the increase of the usage of complex and sophisticated financial tools, the surface become curved. In this sense, both Stochastic Volatility/Jump Diffusion and Different opinions on future trend are not the critical idea for this effect since even before 1987, jump occurs frequently and volatility is still not a constant (eg, volatility cluster).

In this sense, the evolution of volatility surface from flat plane to curved one mainly comes from the hedging side boosted by quants. If someone(some companies) began to use the sophisticated model (eg, Garch model,local volatility models, SV/JD models) instead of BS model to model the options, the volatility surface appear naturally. And then more and more people noticed this effect and then they began to use new tools and then .... It is a feedback loop like electronic amplification circuit.

Within this philosophy, how can we model the volatility surface? The answer is that we need adaptive model can survive as the volatility surface itself is emergent effect due to market feedback of human knowledge. But that can not solve anything. Our model does not have intelligence, they can not evolve dynamically ...

So here, let me just analyze those basic models one by one , and thinking is still going on...

1) local volatility model
The most disadvantage of this model is that the asymptotic surface when maturity time period is large is flat. The reason is simple, when the time period is large, one can go almost all kind of path of different price of stocks, then there is no advantage of usage of local volatility. According to ito33.com, we can not local "everything", interest rate, volatility, correlation can not be parameterized at all, they are dynamically interacted and changed with the evolution of markets.

2) Stochastic volatility model
At first time, it is promising as it is adaptive in some degree. But first, maybe the fatal point, is that option is just the reflection of the human prediction for the market in the future. It is not true that everyone use the SV model or use simply vanilla option to calibration. In this sense, people's prediction can not agree each other toward a option. Secondly, jump happens. One can not mimic large jump with the usage of brownian motion. In mathematic language is that one can compensate the small jump by correcting the model but can not compensate the large jump (That's why we need compensated poisson process to mimic a levy process). It is true that Stochastic volatility can correctly reflect volatility cluster in different economical period but the jump nature determined SV model is not a perfect one . Third, it is not easy to hedge out the risk of uncertainty of volatility.

3) Jump diffusion model
Despite that jump model cannot predict any real jump like the impact of 911, the jump model statistically model the jump in a reasonable sense. And they can correctly mimic the fat tail effect and also fit the distribution well like other sophisticated models.

(to be continued...)

Note03: How can we add our own function to QuantLibAddin

1. Add a new function to an existing category
A. Eg: Our function contains in
QuantLibAddin/qla/****.cpp
QuantLibAddin/qla/****.hpp
B. File generalutils.hpp in the qla directory can be included to pick up additional utility functions.
C. QuantLibAddin is itself a C++ Addin which can be loaded directly to standalone C++ client applications. So it's best to test the new functionality in a standalone program before autogenerating the source for the spreadsheets.
Then Edit the file
QuantLibAddin/Clients/C++/instruments.cpp
adding some example code to demonstrate the use of the new function.
D. Edit QuantLibAddin/srcgen/instruments.xml to provide definition of the new function. See the link for the details of the Edition!

E. Rebuild the srcgen project to generate the source for the Addins
F. Rebuild the Addins.
G. Amend the client files ...
QuantLibAddin/Clients/Excel/instruments.xls QuantLibAddin/Clients/Calc/instruments.sxc QuantLibAddin/Clients/C/instruments.c
... to demonstrate the use of the new function.

2. Add a new category
see link...

Note 02: Seven Easy Steps to debug Excel Add-ins in Visual Studio 2005

Before any step, you should have a program to debug!

Step 1: To build the Debug version, use Visual Studio's Build - Configuration Manager... menu option and select the Debug configuration as Debug mode

Step 2: Select the Project in the Solution Explorer and show its Properties page. Select the "Debugging" node in the left pane, and change the "Command" property to the preferred version of EXCEL.EXE. (You can omit this step if you do not need to change the Excel version from the beginning)

Step 3: Press F5 (or use the Start - Debug menu option) to start Excel.

Step 4: As soon as Excel starts, it will attempt to open your the Debug build of your add-in. (This is because the XLL+ AppWizard set the "Command Arguments" property of the configuration to "$(TargetDir)/Tutorial1.xll" when the project was created.)

Step 5: Use the add-in function in Excel

Step 6: You need to set a break-point in VC 2005 (use the F9 key (or the right-mouse menu) to add a break-point to the add-in function)

Step 7: Let the formula cell calculate. The ways to go to the cell and click on the F2 key then the Enter key. (Alternatively, you can simply change the value of dependent cell and the Excel will automatically recalculate)

After you press Enter, Excel will call the add-in function and the DevStudio debugger will break at the break-point. You can step through the lines of the function (with F10 and F11) or continue (with F5).

Close Excel and return to Visual Studio.

Note 01: Typical Example for using QuantLibAddin C++ client

#include  #include 

using namespace std;
using namespace QuantLib;
using namespace ObjHandler;
using namespace QuantLibAddin;

int main() {
// deal with errors
try {
// setup log file
setLogFile("*****.log");
// direct log messages to stdout
setConsole(1);
// write log message to log file
logMessage("************");

ObjHandler::obj_ptr MyObject(new *****(***));...

// store "MyObject"
storeObject("MyObject",Myobject) ;
//ObjectHandler::Repository::instance().storeObject("MyObject", Myobject);

// write to log file the object's ID
logObject("Myobject");

// retrieve "MyObject" Retrieve the object with the
//given ID and recast it to the desired type
// retrieveObject("MyObJect",Myobject);

return 0;
} catch (const std::exception &e) {
std::ostringstream s;
s << "Error: " << e.what();
ObjectHandler::logMessage(s.str(), 1); return 1;
} catch (...)
{ ObjectHandler::logMessage("Error", 1);
return 1; }

Planning in this week

Get familar with QuantLib, Boost, the way to export functions from QuantLib to excel and the way to debug through Excel to VC8

特殊既得利益集团

“中国特色社会主义是中国共产党人在长期的社会主义建设和改革过程中逐步探索形成的，同民主社会主义在理论和实践上都有着本质的区别。中国特色社会主义始终坚持马克思主义的指导地位，决不搞指导思想多元化，并结合新的实际不断推进马克思主义中国化，用发展着的马克思主义指导新的实践；坚持中国特色社会主义的政治发展道路，坚持中国共产党的领导，决不搞西方的三权分立和多党制；坚持社会主义市场经济体制，坚持社会主义公有制为主体、多种所有制经济共同发展的基本经济制度，坚持按劳分配为主体、多种分配方式并存的分配制度；坚持社会主义先进文化的前进方向，建设社会主义核心价值体系。” （见《人民日报》 ( 2007-05-10 第09版 )）

An introduction to Ox Language( Draft: Mainly copy from others without edition)

Ox is a very easy language if you are familiar with C++. So I believe we can grasp it
in an hour.

A. What is Ox?
Ox is an object-oriented matrix language with a comprehensive mathematical and statistical function library. Matrices can be used directly in expressions, for example to multiply two matrices, or to invert a matrix. The basic syntax elements of Ox are similar to the C, C++ and Java languages (however, knowledge if these languages is not a prerequisite for using Ox). This similarity is most clear in syntax items such as loops, functions, arrays and classes. A major difference is that Ox variables have no explicit type, and that special support for matrices is available.

B. A first look at Ox program
#include // include Ox standard library header

main() // function main is the starting point
{
decl m1, m2; // declare two variables, m1 and m2
//Variables may be declared by using the decl statement, and have no type until the program is actually run.

m1 = unit(3);// assign to m1 a 3 x 3 identity matrix
m1[0][0] = 2; // set top-left element to 2
m2 = <0,0,0;1,1,1>;//m2 is a 2 x 3 matrix, the first
// row consists of zeros, the second of ones
//<0,0,0;1,1,1> is a matrix constant. Elements are listed by row, whereby rows are separated by a semicolon, and elements // within a row by a colon. This value is stored in m2, which is now also of type matrix.
print("two matrices", m1, m2); // print the matrices
}

The program is in ox\samples\myfirst.ox; running this program should produce the following result:


two matrices
2.0000 0.00000 0.00000
0.00000 1.0000 0.00000
0.00000 0.00000 1.0000

0.00000 0.00000 0.00000
1.0000 1.0000 1.0000

C. Variables

Variables are declared using the decl keyword. Unlike C, variables are implicitly typed. This means that variables do not have a type when they are declared, but get a type when values are assigned to them. So a variable can change type during its lifetime. The most important implicit types are int for an integer value, double for a real number, string for a text string and matrix for a matrix (two-dimensional array) of real numbers. The next Ox program illustrates implicit declaration and scope:


#include
main()
{
decl i, d, m, s;
i = 1; // assign integer to i --> i is of type int
d = 1.0; // assign real number to d --> d is double
s = "some text";//assign string to s --> s is string
m = zeros(3,3);/ assign to m a 3 x 3 matrix of zeros
// --> m is of type matrix
print("i=", i, " d=", d, " s=", s, "\nm=", m);
}
This prints (\n is the newline character):


i=1 d=1 s=some text
m=
0.00000 0.00000 0.00000
0.00000 0.00000 0.00000
0.00000 0.00000 0.00000
The scope of a variable refers to the parts of the program which can see the variable. This could be different from its lifetime: a variable can be `alive' but not `seen'. If a variable is declared outside any function, its scope is the remainder of the source file. It is possible to export such variables to other source files, as we shall see shortly.

Variables declared inside a function have scope until the closing brace of the level at which it is declared. The following example illustrates:


#include
decl mX; // external variable
main()
{
decl i = 0; // local variable
{
decl i = 1, j = 0; // new i
mX = ones(3,3);
print("i=", i, " j=", j); // prints: i=1 j=0
} // brace end: local i and j cease to exist
print("\ni=", i); // revert to old i, prints: i=0
}
The variable mX (here we use Hungarian notation), can be seen everywhere in the main function. To make sure that it can never be seen in other source files, prefix it with the word static. It is good programming practice to use static in such cases, because it is very useful to know that it is not used in any other files (we may than rename it, e.g., without any unexpected side effects). An example will be given in myfunc.ox below.


Indexing matrices

Indexing starts at zero, so m[0][0] is the first element of the matrix m. It is easy to select individual elements or a subset of the matrix. Here are some examples:


#include
main()
{
decl m = <0 1 2; 3 4 5; 6 7 8>;
println("m = ", m);
println("element 1,0: ", m[1][0]);
println("second row: ", m[1][]);
println("second column: ", m[][1]);
println("without 1st row/3rd col: ", m[1:][:1]);
println("indexed as a vector ", m[2:3]);
}
Which prints as output:

m =
0.00000 1.0000 2.0000
3.0000 4.0000 5.0000
6.0000 7.0000 8.0000
element 1,0: 3
second row:
3.0000 4.0000 5.0000
second column:
1.0000
4.0000
7.0000
without 1st row/3rd col:
3.0000 4.0000
6.0000 7.0000
indexed as a vector
2.0000
3.0000
These expressions may also be used in assignments, for example:


m[1:][:1] = 10;
m[2:3] = m[6:7];

Functions and function arguments

We have already used various functions from the standard library (such as print, ones and zeros), and written various main functions). Indeed, an Ox program is primarily a collection of functions. It is important to know that all function arguments are passed by value. This means that the function gets a copy which it can change without changing the original. For example:


#include
func(mA)
{
mA = zeros(1,2);
print("ma in func()", mA);
}
main()
{
decl ma;
ma = ones(1,2);
print("ma before func()", ma);
func(ma);
print("ma after func()", ma);
}
which prints:


ma before func()
1.0000 1.0000
ma in func()
0.00000 0.00000
ma after func()
1.0000 1.0000
If the function argument is not changed by the function, it is good programming style to prefix it with the const keyword, as in:


func(const mA)
{
print("ma in func()", mA);
}
Then the compiler can generate much more efficient code, especially for matrices and strings.

Of course it is possible to return changed values from the function. If there is only one return value, this is most simply done by using the return statement:


#include
func(const r, const c)
{
return rann(r, c); // return r x c matrix of random
} // numbers from standard normal
main()
{
print("return value from func():", func(1,2) );
}
Another way is to pass a reference to the variable, rather than the variable itself, as for example in:


#include
func(const pmA)
{
pmA[0] = zeros(1,2);
print("ma in func()", pmA[0]);
}
main()
{
decl ma;
ma = ones(1,2);
print("ma before func()", ma);
func(&ma);
print("ma after func()", ma);
}
which prints:


ma before func()
1.0000 1.0000
ma in func()
0.00000 0.00000
ma after func()
0.00000 0.00000
Now the change to ma is permanent. The argument to the function was the address of ma, and func received that address as a reference. Now we can modify the contents of the reference by assigning a value to pmA[0]. When func has finished, ma has been changed permanently. Note that we gave the argument a const qualification. This was possible because we did not change pmA itself, but what it referred to.


The for and while loops

Since Ox is a matrix language, there is much less need for loop statements than in C or C++. Indeed, because Ox is compiled and then interpreted, there is a speed penalty for using loop statements when they are not necessary.

The for, while and do while loops have the same syntax as in C. The for loop consists of three parts, an initialization part, a termination check, and an incrementation part. The while loops only have a termination check.


#include
main()
{
decl i, d;
for (i = 0; i < 5; ++i)
{
d = i * 0.01;
println(d);
}
}
which prints (println is like print, but ensures that the next output will be starting on a new line):


0
0.01
0.02
0.03
0.04
This could also be written, less elegantly, using while as follows:


#include
main()
{
decl i, d;
i = 0;
while (i < 5)
{
d = i * 0.01;
println(d);
++i;
}
}
It is not uncommon to have more than one loop counter in the for statement, as the following code snippet illustrates:


decl i, j;
for (i = 0, j = 10; i < 5 && j > 0; ++i, --j)
println(i * j);
The && is logical-and, whereas || is logical-or. The ++i statement is called (prefix) incrementation, and means `add one to i'. Similarly, --j subtracts one from j. There is a difference between prefix and postfix incrementation (decrementation). For example, the second line in


i = 3;
j = ++i;
means: add one to i, and assign the result to j, which will get the value 4. But


i = 3;
j = i++;
means: leave the value of i on the stack for assignment, then afterwards increment i. So j will get the value 3. In the incrementation part of the for loop it does not matter whether you use the prefix or postfix form.

D. The if statement

The if statement allows for conditional program flow. In the following example we draw a uniform random number. Such a random number is always between zero and one. The ranu returns a matrix, unless we ask it to generate just one number. Then it returns a double, as is the case here.


decl d = ranu(1,1);
if (d < 0.5)
println("less than 0.5");
else if (d < 0.75)
println("less than 0.75");
else
println("greater than 0.75");
Again, braces are used to group multiple statements together. They should also be used when nesting if statements, to avoid confusion about which else belongs to which if.


decl d1 = ranu(1,1), d2 = ranu(1,1);
if (d1 < 0.5)
{ println("d1 is less than 0.5");
}
else
{ if (d2 < 0.75)
println("d1 >= 0.5 and d2 < 0.75");
else
println("d1 >= 0.5 and d2 <= 0.75");
}
The if part is executed if the expression evaluates to a non-zero value (true). The else part otherwise, i.e. when the expression evaluates to zero (false: either an integer 0, or a double 0.0). Some care is required when using matrices in if statements. A matrix expression is a true statement if all elements are true (non-zero). Even if only one element is zero, the matrix expression is false, so


#include
main()
{
if (ones(2,2)) print("yes");
else print("no");
if (unit(2)) print("yes");
else print("no");
if (zeros(2,2)) print("yes");
else print("no");
}
prints: yesnono.

There are two forms of relational operators. There is < <= > >= == != meaning `less', `less than or equal', `greater', `greater than or equal', `is equal' and `is not equal'. These always produce the integer value 1 (true) or 0 (false). If any of the arguments is a matrix, the result is only true if it is true for each element:


#include
main()
{
if (ones(2,2) == 1) print("yes"); // true for each
else print("no"); // element
if (unit(2) == 1) print("yes");//not true for each
else print("no"); // element
if (zeros(2,2) == 1) print("yes");//not true for each
else print("no"); // element
}
prints: yesnono.

The second form are the dot-relational operators .< .<= .> .>= .== .!= meaning `dot less', `dot less than or equal', `dot greater', `dot greater than or equal', `is dot equal' and `is not dot equal'. If any of the arguments is a matrix, the result is a matrix of zeros and ones, with each element indicating the relevant result.

The any library function returns 1 (true) if any element of the matrix is non-zero, so that yesyesno will be printed by:


#include
main()
{
if (any(ones(2,2))) print("yes");
else print("no");
if (any(unit(2))) print("yes");
else print("no");
if (any(zeros(2,2))) print("yes");
else print("no");
}
To conclude: you can test whether all elements of a matrix m are equal to one (say) by writing: if (m == 1). To test whether any element is equal to one: if (any(m .== 1)). The expression if (m != 1), on the other hand, is only true if none of the elements is equal to one. So, use if (!(m == 1)) to test whether it is true that not all elements are equal to one.


Operations and matrix programming

To a large extent, the same operators are available in Ox as in C or C++. Some of the additional operators are power (^), horizontal concatenation (~), vertical concatenation (|) and the Kronecker product (**). One important distinction is that the operators are also available for matrices, so that, for example, two matrices can be added up directly. For some operators, such as multiplication, there is a distinction between the dot operators (e.g. .* is element by element multiplication and * is matrix multiplication if both arguments are matrices). Not available in Ox are the bitwise operators, instead you need to use the library functions binand and binor.

Because Ox is implicitly typed, the resulting type of the expression will depend on the types of the variables in the expression. When a mixture of types is involved, the result is promoted upwards in the order integer, double, matrix. So in an expression consisting if an integer and a double, the integer will be promoted to a double. An expression of only integers yields an integer. However, there are two important exceptions to this rule:


integer division is done in floating point and yields a double. This is an important difference with C, where integer division is truncated to an integer.
power expressions involving integers which yield a result too large to be expressed as an integer give a double result.

To illustrate, we write the Fahrenheit to Celsius example of Kernighan and Ritchie (1988) in Ox:


#include
const decl LOWER = 0;
const decl UPPER = 100;
const decl STEP = 20;
main()
{
decl fahr;
for (fahr = LOWER; fahr <= UPPER; fahr += STEP)
print("%3d", fahr, " ",
"%6.1f", (5.0/9.0) * (fahr-32), "\n");
}
which prints:


0 -17.8
20 -6.7
40 4.4
60 15.6
80 26.7
100 37.8
In C we have to write 5.0/9.0, because 5/9 evaluates to zero. In Ox both expressions would be evaluated in floating point arithmetic.

In general we get more more efficient code by vectorizing each program as much as possible:


#include
const decl LOWER = 0;
const decl UPPER = 100;
const decl STEP = 20;
main()
{
decl fahr;
fahr = range(LOWER, UPPER, STEP)';
print("%6.1f", fahr ~ (5.0/9.0) * (fahr-32) );
}
As in the first version of the program, we declare three constants which define the Fahrenheit part of the table.

The range() function creates a 1×n matrix with the values LOWER, LOWER+STEP, LOWER + 2STEP, ..., UPPER.

The transpose operator ' changes this into an n×1 matrix.

The conversion to Celsius in the print statement works on the matrix as a whole: multiplication of a matrix by a scalar is equivalent to multiplication by the scalar of each element of the matrix.

The ~ operator concatenates the two column vectors into an n ×2 matrix.

Finally, the print function is different from the printf in C. In Ox each variable to print is simply specified sequentially. It is possible, as done here with "%6.1f", to insert formatting strings for the next variable.

The program prints a table similar to the earlier output:


0.0 -17.8
20.0 -6.7
40.0 4.4
60.0 15.6
80.0 26.7
100.0 37.8

Arrays

The Ox syntax allows for arrays, so you may use, for example, an array of strings (often useful), an array of matrices, or even an array of an array of matrices (etc.). The following program gives an example.


#include
const decl MX_R = 2;
const decl MX_C = 2;
main()
{
decl i, asc, asr, m;
asr = new array[MX_R];
asc = new array[MX_C];
for (i = 0; i < MX_R; ++i)
asr[i] = sprint("row ", i);
for (i = 0; i < MX_C; ++i)
asc[i] = sprint("col ", i);
m = ranu(MX_R, MX_C);
print("%r", asr, "%c", asc, m);
}
which prints


col 0 col 1
row 0 0.020192 0.68617
row 1 0.15174 0.74598
The new operator declares a new object. That could be a class object, as discussed in the next chapter, a matrix, a string, or, as used here, an array. The argument in square brackets is the size of the array. (When creating a matrix in this way, note that a matrix is always two-dimensional, and needs two arguments, as in: m = new matrix[2][2].)

The sprint functions return a string, which is stored in the arrays.

In print(), we use "%r" followed by an array of strings to specify row labels for the subsequent matrix. Columns labels use "%c".


Multiple files: using #include and #import

The source code of larger projects will often be spread over several source files. Usually the .ox file containing the main function is only a few tens of lines. We have already seen that information about other source files is passed on through included header files. However, to run the entire program, the code of those files needs to be linked together as well. Ox offers various ways of doing this. As an example, consider a mini-project consisting of two files: a source code file and a header file. The third file will contain the main function.

Mini project [samples/myfunc.ox]

#include
static decl s_iCalls = 0; // counter, initialize to 0
MyFunction(const ma)
{
++s_iCalls; // increment calls counter
println("MyFunction has been called ", s_iCalls,
" times and prints:", ma);
}
[samples/myfunc.h]

MyFunction(const ma);
The header file myfunc.h declares the MyFunction function, so that it can be used in other Ox files. Note that the declaration ends in a semicolon. The source code file contains the definition of the function, which is the actual code of the function. The header of the definition does not end in a semicolon, but is followed by the opening brace of the body of the function. The s_iCalls variable is declared outside any function, making it an external variable. Here we also use the static type specifier, which restricts the scope of the variable to the myfunc.ox file: s_iCalls is invisible anywhere else (and other files may contain their own s_iCalls variable). Variables declared inside a block of curly braces have a more limited lifetime. Their life starts when they are declared, and finishes at the closing brace (matching the brace level of declaration).

It is also possible to share variables between various source files, although there can be only one declaration (physical allocation) of the shared variable. The following modifications would do that for the myfunc.ox program:
(1) delete the static keyword from the declaration,
(2) add to myfunc.h the line (renaming s_iCalls to g_iCalls):


extern decl g_iCalls;
Any code which includes myfunc.h can now reference or change the g_iCalls variable.


Including the code into the main file

The first way of combining the mini project with the main function is to #include the actual code. In that case the myfunc.h header file is not needed:

[samples/mymaina.ox]

#include
#include "myfunc.ox"
main()
{
MyFunction("one");
}
The result will be just one code file, and mymaina.ox can be run as oxl mymaina.


Importing the code into the main file

The drawback of the previous method of including source code using #include, is that it can only be done once. That is not a problem in this short program, but is difficult to ensure if a library is used at many points in a large project. The #import command solves this problem.

[samples/mymainc.ox]

#include
#import "myfunc"
main()
{
MyFunction("one");
}
Again, mymainc.ox can be run as oxl mymainc. There is no extension in the argument to #import. The effect is as an #include "myfunc.h" statement followed by marking myfunc.ox for linking.
[ #import will actually try to find the .oxo file first. If that is not found, it will search for the .ox file. If neither is found, the program cannot run. More detail is in import of modules. ]
The actual linking only happens when the file is run, and the same #import statement may occur multiple times (as well as in compiled files). So even when the same file is imported many times, it will only be linked once.


Importing Ox packages

If myfunc.ox would require the maximization package, it could have at the top:


#include
#import
Partial paths can be used. Searching is relative to the OX4PATH environment variable. For example, if the Arfima package is in its default location of ox/packages/arfima, we would use:


#import
The distinction between angular brackets and double quotes in the include and import statements is discussed in file inclusion. Roughly, the <> form should be used for files which are part of the Ox system, and the double quotes for your own files, which will not be in the Ox tree.


Separate compilation

Ox source code files can be compiled into Ox object files. These files have the .oxo extension, and are binary. The format is identical across operating systems, but since they are binary, transfer from one platform to another has to be done in binary mode.

To compile myfunc.ox into an Ox object file use the -c switch:


oxl -c myfunc
This creates myfunc.oxo (the .oxo extension is automatically appended). Remember that myfunc.oxo must be recreated every time myfunc.ox changes.

Now, when rerunning mymainc.ox, it will automatically use the .oxo instead of the .ox file.

Compiled Ox files can be useful for very large files (although even then compilation will be very fast), or if you do not wish to distribute the source files. They are inconvenient when the code is still under development.


Object-oriented programming

Object-oriented programming involves the grouping together of functions and variables in convenient building blocks. These blocks can then be used directly, or as starting point for a more specialized implementation. A major advantage of object-oriented programming is that it avoids the use of global variables, thus making the code more re-entrant: several instances will not conflict wiith each other.

The object-oriented features in Ox are not as sophisticated as in some low-level languages. However, this avoids the complexity of a language such as C++, while still providing most of the benefits.

Ox allows you to completely ignore the object-oriented features. However, you will then not be able to use the preprogrammed classes for data management and simulation. It is especially in the latter task that we found a considerable reduction in the required programming effort after writing the base class.

The class is the main vehicle for object-oriented programming. A class is nothing more than a group of variables (the data) and functions (the actions) packaged together. This makes it a supercharged struct (or record in Pascal terminology). Inheritance allows for a new class to add data and functions to the base class, or even redefine functionality of the base class.

In Ox, the default is that all data members of the class are protected (only visible to class members), and all function members are public. Like C++, Ox has the virtual keyword to define functions which can be replaced by the derived class. Classes are used by dynamically creating objects of that class. No static objects exist in Ox. When an object is created, the constructor function is called, when the object is deleted, the destructor function is called. More information on object-oriented programming is given in the syntax chapter.


Style and Hungarian notation

The readability and maintainability of a program is considerably enhanced when using a consistent style and notation, together with proper indentation and documentation. Style is a personal matter; this section describes the one I have adopted.

In my code, I always indent by one tab (four spaces) at the next level of control (i.e. after each opening brace), jumping back on the closing brace.


Table tut.1: Hungarian notation prefixes

prefix type example
i integer iX
c count of cX
b boolean (f is also used) bX
fl integer flag flX
d double dX
m matrix mX
v vector (1 ×n or n ×1 matrix) vX
s string sX
as array of strings asX
am array of matrices amX
a reference in function argument amX
m_ class member variable m_mX
s_ static external variable (file scope) s_mX
g_ external variable with global scope g_mX
fn function reference fnX
I have found Hungarian notation especially useful (see e.g. Petzold, 1992, Ch. 1). Hungarian notation involves the decoration of variable names. There are two elements to Hungarian notation: prefixing of variable names to indicate type (Table Table 1), and using case to indicate scope (Table Table 2, remember that Ox is case sensitive).


Table tut.1: Hungarian notation, case sensitivity

function all lowercase
function (exported) first letter uppercase
static external variable type in lowercase, next letter uppercase
(perhaps prefixed with s_)
exported external variable as above, but prefixed with g_
function argument type in lowercase, next letter uppercase
local variables all lowercase
constants all uppercase
As an example consider:


#include
const decl MX_R = 2; /* a constant */
decl g_mX; /* exported matrix */
static decl s_iCount; /* static external variable */
static func1(const pdX)/* argument is pointer to double */
{
}
/* exported function */
Func2(const mX, const asX, const cT, const cX)
{
decl i, m;
}
Func2 expects a cT × cX matrix, and corresponding array of cX variable names. The c prefix is used for the number of elements in a matrix or string. Note however, that it is not necessary in Ox to pass dimensions separately. You can ask mX and asX what dimensions they have:


Func2(const mX, const asX)
{
decl i, m, ct, cx;
cx = columns(mX);
ct = rows(mX);
if (cx != sizeof(asX))
print("error: dimensions don't match");
}

Optimizing for speed

Ox is very fast: current benchmarks suggest that it is faster than most (if not all) other commonly used matrix language interpreters. A program can never be fast enough though, and here are some tips to achieve even higher speed:


Use matrices as much as you can, avoiding loops and matrix indexing.
Use the const argument qualifier when an argument is not changed in a function: this allows for more efficient function calling.
Use built-in functions where possible.
When optimizing a program with loops, it usually only pays to optimize the inner most loop. One option is to move loop invariants to a variable outside the loop.
Avoid using `hat' matrices, i.e. avoid using outer products over large dimensions when not necessary.
Note that matrices are stored by row (the C and C++ default, but transposed from the Fortran default), so it could sometimes be faster to transpose matrices (i.e. have data variables in rows instead of columns).
If necessary, you can link in C or Fortran code.


OxGauss
Ox has the capability of running a wide range of Gauss (GAUSS is a trademark of Aptech Systems, Inc., Maple Valley, WA, USA) programs. Gauss code can be called from Ox programs, or run on its own. The formal syntax of OxGauss is described in the Ox Appendices, which also lists some of the limitations of OxGauss, and gives a function summary. The remainder of this chapter gives some examples on its use. More information on OxGauss can be found:

in the Ox Appendices,
on the web: Running GAUSS programs Under Ox3 by Philip A. Viton.
Running OxGauss programs from the command line
As an example we consider a small project, consisting of a code file that contains a procedure and an external variable, together with a code file that includes the former and calls the function. We shall use the .src or .oxgauss extension for the OxGauss programs.

[samples/oxgauss/gaussfunc.src]


declare matrix _g_base = 1;

proc(0)=gaussfunc(a,b);
"calling gaussfunc";
retp(a+_g_base*eye(b));
endp;
[samples/oxgauss/gausscall.src]


#include gaussfunc.src;

_g_base = 20;
z = gaussfunc(10,2);
"result from gaussfunc" z;
To run this program on the command line, enter


oxl -g gausscall.src
Which produces the output:


Ox version 4.00 (Windows/U) (C) J.A. Doornik, 1994-2006
calling gaussfunc
result from gaussfunc
30.000000 10.000000
10.000000 30.000000
If there are problems at this stage, we suggest to start by reading the first chapter of the Introduction to Ox.

Running OxGauss programs from OxMetrics
Using Ox Professional, the OxGauss program can be loaded into OxMetrics. The syntax highlighting makes understanding the program easier. Click on Run (the running person) to execute the program. This runs the program using the OxGauss application, with the output in a window entitled `OxGauss Session'. OxMetrics will treat the file as an OxGauss file if it has the .src, .g or .oxgauss extension. If not, the file can still be run by launching OxGauss from the OxMetrics workspace window.


Calling OxGauss from Ox
The main objective of creating OxGauss was to allow Gauss code to be called from Ox. This helps in the transition to Ox, and increases the amount of code that is available to users of Ox. The main point to note is that the OxGauss code lives inside the gauss namespace. In this way, the Ox and OxGauss code can never conflict. Returning to the earlier example, the first requirement is to make an Ox header file for gaussfunc.src. This must declare the external variables and procedures explicitly in the gauss namespace:

[samples/oxgauss/gaussfunc.h]


namespace gauss
{
extern decl _g_base;
gaussfunc(const a, const b);
}
Next, the OxGauss code must be imported into the Ox program. The #import command has been extended to recognize OxGauss imports by prefixing the file name with gauss::, as in the following program:

[samples/oxgauss/gausscall.ox]


#include
#import "gauss::gaussfunc"
main()
{
gauss::_g_base = 20;
decl z = gauss::gaussfunc(10,2);
println("result from gaussfunc", z);
}
When the OxGauss functions or variables are accessed, they must also be prefixed with the namespace identifier gauss::. The output is:


calling gaussfunc
result from gaussfunc
30.000 10.000
10.000 30.000

A good Random number generator (C++ code)

I recommend the following code written by Richard J. Wagner.

// MersenneTwister.h
// Mersenne Twister random number generator -- a C++ class MTRand
// Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
// Richard J. Wagner v1.0 15 May 2003 rjwagner@writeme.com

// The Mersenne Twister is an algorithm for generating random numbers. It
// was designed with consideration of the flaws in various other generators.
// The period, 2^19937-1, and the order of equidistribution, 623 dimensions,
// are far greater. The generator is also fast; it avoids multiplication and
// division, and it benefits from caches and pipelines. For more information
// see the inventors' web page at http://www.math.keio.ac.jp/~matumoto/emt.html

// Reference
// M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-Dimensionally
// Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions on
// Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3-30.

// Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
// Copyright (C) 2000 - 2003, Richard J. Wagner
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions
// are met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. The names of its contributors may not be used to endorse or promote
// products derived from this software without specific prior written
// permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

// The original code included the following notice:
//
// When you use this, send an email to: matumoto@math.keio.ac.jp
// with an appropriate reference to your work.
//
// It would be nice to CC: rjwagner@writeme.com and Cokus@math.washington.edu
// when you write.


// Not thread safe (unless auto-initialization is avoided and each thread has
// its own MTRand object)

#include
#include
#include
#include
#include

class MTRand {
// Data
public:
typedef unsigned long uint32; // unsigned integer type, at least 32 bits

enum { N = 624 }; // length of state vector
enum { SAVE = N + 1 }; // length of array for save()

protected:
enum { M = 397 }; // period parameter

uint32 state[N]; // internal state
uint32 *pNext; // next value to get from state
int left; // number of values left before reload needed


//Methods
public:
MTRand( const uint32& oneSeed ); // initialize with a simple uint32
MTRand( uint32 *const bigSeed, uint32 const seedLength = N ); // or an array
MTRand(); // auto-initialize with /dev/urandom or time() and clock()

// Do NOT use for CRYPTOGRAPHY without securely hashing several returned
// values together, otherwise the generator state can be learned after
// reading 624 consecutive values.

// Access to 32-bit random numbers
double rand(); // real number in [0,1]
double rand( const double& n ); // real number in [0,n]
double randExc(); // real number in [0,1)
double randExc( const double& n ); // real number in [0,n)
double randDblExc(); // real number in (0,1)
double randDblExc( const double& n ); // real number in (0,n)
uint32 randInt(); // integer in [0,2^32-1]
uint32 randInt( const uint32& n ); // integer in [0,n] for n < 2^32
double operator()() { return rand(); } // same as rand()

// Access to 53-bit random numbers (capacity of IEEE double precision)
double rand53(); // real number in [0,1)

// Access to nonuniform random number distributions
double randNorm( const double& mean = 0.0, const double& variance = 0.0 );

// Re-seeding functions with same behavior as initializers
void seed( const uint32 oneSeed );
void seed( uint32 *const bigSeed, const uint32 seedLength = N );
void seed();

// Saving and loading generator state
void save( uint32* saveArray ) const; // to array of size SAVE
void load( uint32 *const loadArray ); // from such array
friend std::ostream& operator<<( std::ostream& os, const MTRand& mtrand );
friend std::istream& operator>>( std::istream& is, MTRand& mtrand );

protected:
void initialize( const uint32 oneSeed );
void reload();
uint32 hiBit( const uint32& u ) const { return u & 0x80000000UL; }
uint32 loBit( const uint32& u ) const { return u & 0x00000001UL; }
uint32 loBits( const uint32& u ) const { return u & 0x7fffffffUL; }
uint32 mixBits( const uint32& u, const uint32& v ) const
{ return hiBit(u) | loBits(v); }
uint32 twist( const uint32& m, const uint32& s0, const uint32& s1 ) const
{ return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL); }
static uint32 hash( time_t t, clock_t c );
};

inline MTRand::MTRand( const uint32& oneSeed )
{ seed(oneSeed); }

inline MTRand::MTRand( uint32 *const bigSeed, const uint32 seedLength )
{ seed(bigSeed,seedLength); }

inline MTRand::MTRand()
{ seed(); }

inline double MTRand::rand()
{ return double(randInt()) * (1.0/4294967295.0); }

inline double MTRand::rand( const double& n )
{ return rand() * n; }

inline double MTRand::randExc()
{ return double(randInt()) * (1.0/4294967296.0); }

inline double MTRand::randExc( const double& n )
{ return randExc() * n; }

inline double MTRand::randDblExc()
{ return ( double(randInt()) + 0.5 ) * (1.0/4294967296.0); }

inline double MTRand::randDblExc( const double& n )
{ return randDblExc() * n; }

inline double MTRand::rand53()
{
uint32 a = randInt() >> 5, b = randInt() >> 6;
return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0); // by Isaku Wada
}

inline double MTRand::randNorm( const double& mean, const double& variance )
{
// Return a real number from a normal (Gaussian) distribution with given
// mean and variance by Box-Muller method
double r = sqrt( -2.0 * log( 1.0-randDblExc()) ) * variance;
double phi = 2.0 * 3.14159265358979323846264338328 * randExc();
return mean + r * cos(phi);
}

inline MTRand::uint32 MTRand::randInt()
{
// Pull a 32-bit integer from the generator state
// Every other access function simply transforms the numbers extracted here

if( left == 0 ) reload();
--left;

register uint32 s1;
s1 = *pNext++;
s1 ^= (s1 >> 11);
s1 ^= (s1 << 7) & 0x9d2c5680UL;
s1 ^= (s1 << 15) & 0xefc60000UL;
return ( s1 ^ (s1 >> 18) );
}

inline MTRand::uint32 MTRand::randInt( const uint32& n )
{
// Find which bits are used in n
// Optimized by Magnus Jonsson (magnus@smartelectronix.com)
uint32 used = n;
used |= used >> 1;
used |= used >> 2;
used |= used >> 4;
used |= used >> 8;
used |= used >> 16;

// Draw numbers until one is found in [0,n]
uint32 i;
do
i = randInt() & used; // toss unused bits to shorten search
while( i > n );
return i;
}


inline void MTRand::seed( const uint32 oneSeed )
{
// Seed the generator with a simple uint32
initialize(oneSeed);
reload();
}


inline void MTRand::seed( uint32 *const bigSeed, const uint32 seedLength )
{
// Seed the generator with an array of uint32's
// There are 2^19937-1 possible initial states. This function allows
// all of those to be accessed by providing at least 19937 bits (with a
// default seed length of N = 624 uint32's). Any bits above the lower 32
// in each element are discarded.
// Just call seed() if you want to get array from /dev/urandom
initialize(19650218UL);
register int i = 1;
register uint32 j = 0;
register int k = ( N > seedLength ? N : seedLength );
for( ; k; --k )
{
state[i] =
state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL );
state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
state[i] &= 0xffffffffUL;
++i; ++j;
if( i >= N ) { state[0] = state[N-1]; i = 1; }
if( j >= seedLength ) j = 0;
}
for( k = N - 1; k; --k )
{
state[i] =
state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL );
state[i] -= i;
state[i] &= 0xffffffffUL;
++i;
if( i >= N ) { state[0] = state[N-1]; i = 1; }
}
state[0] = 0x80000000UL; // MSB is 1, assuring non-zero initial array
reload();
}


inline void MTRand::seed()
{
// Seed the generator with an array from /dev/urandom if available
// Otherwise use a hash of time() and clock() values

// First try getting an array from /dev/urandom
/* FILE* urandom = fopen( "/dev/urandom", "rb" );
if( urandom )
{
uint32 bigSeed[N];
register uint32 *s = bigSeed;
register int i = N;
register bool success = true;
while( success && i-- )
success = fread( s++, sizeof(uint32), 1, urandom );
fclose(urandom);
if( success ) { seed( bigSeed, N ); return; }
} */

// Was not successful, so use time() and clock() instead
seed( hash( time(NULL), clock() ) );
}


inline void MTRand::initialize( const uint32 seed )
{
// Initialize generator state with seed
// See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
// In previous versions, most significant bits (MSBs) of the seed affect
// only MSBs of the state array. Modified 9 Jan 2002 by Makoto Matsumoto.
register uint32 *s = state;
register uint32 *r = state;
register int i = 1;
*s++ = seed & 0xffffffffUL;
for( ; i < N; ++i )
{
*s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL;
r++;
}
}


inline void MTRand::reload()
{
// Generate N new values in state
// Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
register uint32 *p = state;
register int i;
for( i = N - M; i--; ++p )
*p = twist( p[M], p[0], p[1] );
for( i = M; --i; ++p )
*p = twist( p[M-N], p[0], p[1] );
*p = twist( p[M-N], p[0], state[0] );

left = N, pNext = state;
}


inline MTRand::uint32 MTRand::hash( time_t t, clock_t c )
{
// Get a uint32 from t and c
// Better than uint32(x) in case x is floating point in [0,1]
// Based on code by Lawrence Kirby (fred@genesis.demon.co.uk)

static uint32 differ = 0; // guarantee time-based seeds will change

uint32 h1 = 0;
unsigned char *p = (unsigned char *) &t;
for( size_t i = 0; i < sizeof(t); ++i )
{
h1 *= UCHAR_MAX + 2U;
h1 += p[i];
}
uint32 h2 = 0;
p = (unsigned char *) &c;
for( size_t j = 0; j < sizeof(c); ++j )
{
h2 *= UCHAR_MAX + 2U;
h2 += p[j];
}
return ( h1 + differ++ ) ^ h2;
}


inline void MTRand::save( uint32* saveArray ) const
{
register uint32 *sa = saveArray;
register const uint32 *s = state;
register int i = N;
for( ; i--; *sa++ = *s++ ) {}
*sa = left;
}


inline void MTRand::load( uint32 *const loadArray )
{
register uint32 *s = state;
register uint32 *la = loadArray;
register int i = N;
for( ; i--; *s++ = *la++ ) {}
left = *la;
pNext = &state[N-left];
}


inline std::ostream& operator<<( std::ostream& os, const MTRand& mtrand )
{
register const MTRand::uint32 *s = mtrand.state;
register int i = mtrand.N;
for( ; i--; os << *s++ << "\t" ) {}
return os << mtrand.left;
}


inline std::istream& operator>>( std::istream& is, MTRand& mtrand )
{
register MTRand::uint32 *s = mtrand.state;
register int i = mtrand.N;
for( ; i--; is >> *s++ ) {}
is >> mtrand.left;
mtrand.pNext = &mtrand.state[mtrand.N-mtrand.left];
return is;
}



// Change log:
//
// v0.1 - First release on 15 May 2000
// - Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
// - Translated from C to C++
// - Made completely ANSI compliant
// - Designed convenient interface for initialization, seeding, and
// obtaining numbers in default or user-defined ranges
// - Added automatic seeding from /dev/urandom or time() and clock()
// - Provided functions for saving and loading generator state
//
// v0.2 - Fixed bug which reloaded generator one step too late
//
// v0.3 - Switched to clearer, faster reload() code from Matthew Bellew
//
// v0.4 - Removed trailing newline in saved generator format to be consistent
// with output format of built-in types
//
// v0.5 - Improved portability by replacing static const int's with enum's and
// clarifying return values in seed(); suggested by Eric Heimburg
// - Removed MAXINT constant; use 0xffffffffUL instead
//
// v0.6 - Eliminated seed overflow when uint32 is larger than 32 bits
// - Changed integer [0,n] generator to give better uniformity
//
// v0.7 - Fixed operator precedence ambiguity in reload()
// - Added access for real numbers in (0,1) and (0,n)
//
// v0.8 - Included time.h header to properly support time_t and clock_t
//
// v1.0 - Revised seeding to match 26 Jan 2002 update of Nishimura and Matsumoto
// - Allowed for seeding with arrays of any length
// - Added access for real numbers in [0,1) with 53-bit resolution
// - Added access for real numbers from normal (Gaussian) distributions
// - Increased overall speed by optimizing twist()
// - Doubled speed of integer [0,n] generation
// - Fixed out-of-range number generation on 64-bit machines
// - Improved portability by substituting literal constants for long enum's
// - Changed license from GNU LGPL to BSD

关于对中国未来的一点思考（一）

近日读到中国人民大学前副校长谢韬在炎黄春秋2007年第二期发表《民主社会主义模式与中国前途》一文，勾起了自己想写点什么东西的想法，本来这个blog是作为我攻读博士学位研究的一个日志，不想写其他东西，但是还是忍不住动手写点思考，毕竟作为炎黄子孙，对中国未来，不仅切身相关，也关系后代子孙。 如果国人能人人深思未来，择其最优行之，我华夏如何不能收复台湾？如何不能和平崛起？

来美国三年，和这边老美教授，学生们聊天时也时有谈起共产主义，其实也很笑人，几乎和我所有聊天的老美（当然也包括来自世界各地的留学生）都承认共产主义是一个理想的完美社会形态，毕竟按共产主义的目标来说，人人平等，自由的理念不但在政治领域，也在经济领域完全实现。 这种社会几近于天堂！ 但是，正如谢先生在《民主社会主义模式与中国前途》一文自问：“德国人是不是应该比我们更懂得马克思，俄国人是不是应该比我们更懂得列宁，就像我们比外国人更懂得孔夫子一样？”。 这么好的东西为什么外国人不要？ 或者要了，为什么最后也都以失败而告终？

抱读马列主义的老学究这个时候可能跳出来说，那是国际 环境的影响，以美国为首的敌对势力的反扑的结果， 不见得能证明那种制度谁强谁弱，谁优谁劣。 的确，胜败乃兵家常事，前苏联为首的共产主义阵线的瓦解是多方面因素的结果，我个人不想在这个问题上多纠缠，但是如果回顾历史，有讽刺意味的是，瓦解竟然是没有经过战争，而是被内部攻破，从某种意义上说，这种瓦解可以看着在当时特定国际环境，条件下，历史必然的选择的结果了，用愤青的话来说，正所谓得道多助，失道寡助。我相信我们党在这个问题上绝对是下了大力气，大功夫研究和分析的，分析成果也肯定变成了92年后邓公进一步改革的理论基础之一。

但是，改革现在完结了吗？ 我们可以一劳永逸的抱着现在的制度吗？ 其实，任何一个还在思考的中国人，只要他还有良心，还愿意去考虑别人，而不是单以自己最大利益化的那些垃圾，随便走到民间，听听人民大众的声音，也知道这种现状必然不会持久。 这也是现在 ”左派“希望回复公有制， 老毛那种虽然也许一穷二白，但是”真正人民做主人“的时代， ”右派“ 们呼吁进行政治改革，甚至如谢先生所说那样，把民主社会主义模式当做镜子，作为参照。

我在这里，不想对上面两派的看法说太多， 但是个人认为要回复老毛时代那种体制，估计难度相当于封建复辟，不是说不能成功，毕竟这种事情有先例，袁世凯也成功了那么几天，但要回到那种类似于以精神领袖（当然不一定特指某个人）为核心的神权社会， 说白了，和邪教复辟无疑，除非经过大规模再洗脑，现在的人民是万万不会接受的，即使 那位”左派“以超常手段成功，也不过是水中花，镜中月，落得个万世骂名而已。 而右派的建议，很多貌似比较可行，就拿美国的民主体制，实现了这么多年，很多国家也有类似的建制，虽然在俺们学究眼中，可以说充满了无穷的弊端，但如果能够完完整整搬到中国，我敢说不出十年，没有人会再说“这种制度在中国注定失败”的低级，毫无头脑的话。 但问题是他们的制度一定是最适合（注意，是最）中国人民的吗？ 为什么我们要选择这种制度？ 更大的问题当然还是在实践方面，这种制度能够完完整整搬到中国？ 行不通嘛！ 现在民间有句笑话“不改革亡国，改革亡党”的笑话正是上面矛盾的真实写照。

谢先生的这篇文章，其实给了两边台阶下的出路，以修正原来理解的马克思理论的基础上，以别人成功了的民主社会主义模式为蓝本，即不走资本主义道路，也维持了走中国特色社会主义道路的党的地位。 当然即使这样，阻力也是重重。 随便google一下，就有好多评论，虽然网民多半叫好，但也有人用“骗子”的大帽子扣了过来， 比如 关于上海，浙江召开“关于《民主社会主义模式与中国前途》问题座谈会” 的报道就提到，“座谈会上同志们满怀义愤，畅所欲言。”谢先生随便一个具有学术讨论性质的文章，就让我们的同志们满怀义愤，集会迎头痛击“民主社会主义”反动思潮，还准备深入进行清算。 用网络上的话来说，这帮人也真算得上是极品了！

Quasi Monte Carlo Integration

According to the theory, the error of quasi Monte
Carlo integration is about O(N^-1), contrasting
to O(N^-1/2) for ordinary Monte Carlo method.

It is noted that O(N^-1) has huge speed advantage
comparing to Monte Carlo method.

The following is the fortran 90 code for quasi Monte
Carlo method with Halton sequence in 4D. Notice
here we do not use parallel version for Halton
sequence.

MODULE GLOBE

REAL*8,PARAMETER::LAMBDA=3.0d0
REAL*8,PARAMETER::GAMMA=1.0d0
REAL*8,PARAMETER::PI=3.141592653589793238462643383279592884197169399375d0
REAL*8,PARAMETER::PI2=PI*2.0d0

ENDMODULE GLOBE



Program quasimontecarlo
use GLOBE
external F
include 'mpif.h'
REAL*8:: a,b,c(1:4),func,kx,ky,qx,qy
INTEGER*8 ,DIMENSION(1:4) :: N1,NRA1, NRB1
REAL*8,DIMENSION(1:4) :: RS1
REAL*8::RESULTI,x,y,z,G,REALvalue
integer:: myid,numprocs,ierr1,timer0,timer
INTEGER*8::ISTEP,i,j
integer::L
call MPI_INIT( ierr1 )
call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr1 )
call MPI_COMM_SIZE(MPI_COMM_WORLD,numprocs,ierr1)
print*,'Lambda=',Lambda, ' Gamma=',Gamma
L=myid+1
timer0=mpi_wtime()
NRA1(:)=0
NRB1(:)=1
N1(1)=2
N1(2)=3
N1(3)=5
N1(4)=7
RESULTI=0.0d0
ISTEP=50000000000
! print*,'test 2'
DO j=1,20
CALL HALTON( N1,NRA1, NRB1,RS1)
END DO
! print *, 'test 1'

DO j=1,ISTEP
CALL HALTON( N1,NRA1, NRB1,RS1)
kx=RS1(1)*PI2
ky=RS1(2)*PI2
qx=RS1(3)*PI2
qy=RS1(4)*PI2
! print*,x,y,z
RESULTI=RESULTI+F(kx,ky,qx,qy,L)
END DO
RESULTI= RESULTI/dble(ISTEP)
timer=mpi_wtime()
print*, myid+1, RESULTI
timer=mpi_wtime()
print*, myid+1, RESULTI

! print*, 'FOR single ISTEP=',ISTEP, (timer-timer0), 'second with result ',RESULTI

call mpi_finalize(ierr1)

END PROGRAM quasimontecarlo




FUNCTION F(kx,ky,qx,qy,L) RESULT(F_RESULT)
USE GLOBE
Integer:: L
REAL*8::x,y,z,F_RESULT,alpha,f1,f2,kx,ky,qx,qy,g1,g2,o1,o2,mk
REAL*8:: fejer1,fejer2,dl ,dc

! p1=(dsin(x)**2)*(dsin(y)**2)*(dsin(z)**2)
! p2= (dcos(x)**2)*(dcos(y)**2)*(dcos(z)**2)
! The function ATAN2(Y,X) = arctan(y,x) will return a value in (-pi, pi].
! If Y is positive the result will be positive. If Y is zero the result will
! be zero if X is positive, and pi if X is negative. If Y is negative the
! result will be negative. If X is zero the result will be plus or minus pi/2.
! Both X and Y are not permitted to be zero simultaneously. The purpose of
! the function is to avoid division by zero.
dl=dble(L)
f1=LAMBDA-(dcos(kx)+dcos(ky))
f2=LAMBDA-(dcos(kx+qx)+dcos(ky+qy))
g1=GAMMA* (dsin(kx)+dsin(ky))
g2=GAMMA* (dsin(kx+qx)+dsin(ky+qy))
o1=Dsqrt(f1*f1+g1*g1)
o2=Dsqrt(f2*f2+g2*g2)
If (( o1 .eq. 0.0d0 ) .OR. (o2 .eq. 0.0d0)) THEN
mk=1.0d0
else
mk=1.0d0- (f1*f2+g1*g2)/(o1*o2)
end if
dc=Dcos(qx)
If ( dc .eq. 1.0d0) THEN
fejer1=dl*dl
else
fejer1=(Dcos(qx*dl)-1.0d0)/(dc-1.0d0)
end if
dc=Dcos(qy)
If ( dc .eq. 1.0d0) THEN
fejer2=dl*dl
else
fejer2=(Dcos(qy*dl)-1.0d0)/(dc-1.0d0)
end if


F_RESULT=fejer1*fejer2*mk
F_RESULT=fejer1*fejer2*mk

ENDFUNCTION F




!----------------Halton sequence-----------------------
! The initialization is:
!NRA1 = 0
!NRB1 = 1
!N1 = p, the base one wishes to use. the pth prime
!DIM is the dimension
SUBROUTINE HALTON( N1,NRA1, NRB1,RS1)
INTEGER :: i
INTEGER*8 :: NXA1,NXB1
INTEGER*8,DIMENSION(1:4) :: N1,NRA1, NRB1
REAL*8,DIMENSION(1:4) :: RS1
DO i=1,4
NXA1 = NRB1(i) - NRA1(i)
! Special case:
IF (NXA1 .EQ. 1) THEN
NRA1(i) = 1
NRB1(i) = NRB1(i)*N1(i)
GOTO 19
END IF
! General case:
NXB1 = NRB1(i)/N1(i)
10 IF (NXA1 .LE. NXB1) THEN
NXB1 = NXB1/N1(i)
GOTO 10
ELSE
NRA1(i) = (1 + N1(i))*NXB1- NXA1
END IF
19 RS1(i) = dble(NRA1(i))/dble(NRB1(i))
END DO
RETURN
END SUBROUTINE HALTON
!--------------------------------------------

Can we make a bet on stock price ?

In http://www.wilmott.com/messageview.cfm?catid=8&threadid=48034, lord12 asked a simple question that r = 0.1 u = 1.2 d=.9. John thinks the stock will go up and makes a bet with Jerry. If S1>S0, John will $100 from Jerry. If S0 is smaller than S1, John will give $100 to Jerry. Who has the advantage?

My idea is the following:
I think they can hedge out the risk which depends on the situation of option at that time. The clue is :


1) They can not hedge out risk if only hold/sell the stock...

In the side of Jerry:

Jerry can buy X stock(Let us suppose the price is 1) at time 0, if stock is up, he will get X*1.2-100 ,if stock is down, he will get X*0.9+100. So, he need X*1.2-100>(1+r)*x=1.1X and X*0.9+100>1.1X to hedge out any risk, 0.1X>100 and 100>0.2X hold at the same time. it is impossible .

In the side of John,

he need X*1.2+100>1.1X and X*0.9-100>1.1X So X>-1000 and -100>0.2X (X<-500). We get -500<-1000! it is impossible if John buy stock.

2) But it seems promising if they can use call/put option to hedge out the risk. Eg: John can sell the call option at time 0 and the strike is K. the call option price per stock is then (1.2-K)*2/3. Suppose he sell X call option to his client, the total money at time 1 is if stock is up, he will get (1+r)*X*(1.2-K)*2/3-(1.2-K)*X+100. So if 100>(1.2-K)*0.27X, he will win 100- (1.2-K)*0.27X bucks.
Notice we need X<100/((1.2-k)*0.27) 100=""> 0.73*(1.2-K)X+KX-0.9X-100, we need 0.23K*X-0.02X-100>0 to get profit, so he needs X>100/(0.23K-0.02) .

In this way, we know if 100/((1.2-K)*0.27)>100/(0.23K-0.02)=> 0.23K-0.02>0.32-0.27K=> 0.5K>0.34 => K>0.68=> Always hold if we know u=1.2 d=0.9.

In this example, if John is a banker, he can bet with Jerry as he can hedge out the risk completely!

Research Note 4/17/07

To calculate the Tr(I-M) which is a lower bound for the block entropy in bipartite fermion systems in any kind of shape, I find that we need to modify the ordinary Fejer kernel. I will keep thinking in this problem.