diff --git a/QuantLib.vcxproj b/QuantLib.vcxproj
index 5fd1389b66c..589e1f39327 100644
--- a/QuantLib.vcxproj
+++ b/QuantLib.vcxproj
@@ -142,7 +142,7 @@
       <InlineFunctionExpansion>AnySuitable</InlineFunctionExpansion>
       <IntrinsicFunctions>false</IntrinsicFunctions>
       <FavorSizeOrSpeed>Speed</FavorSizeOrSpeed>
-      <AdditionalIncludeDirectories>.;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
+      <AdditionalIncludeDirectories>C:\Program Files\boost\boost_1_61_0;.;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
       <PreprocessorDefinitions>NDEBUG;WIN32;_LIB;_SCL_SECURE_NO_DEPRECATE;_CRT_SECURE_NO_DEPRECATE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
       <StringPooling>true</StringPooling>
       <RuntimeLibrary>MultiThreadedDLL</RuntimeLibrary>
@@ -490,6 +490,7 @@
     <ClInclude Include="ql\experimental\finitedifferences\bsmrndcalculator.hpp" />
     <ClInclude Include="ql\experimental\finitedifferences\fdmhestongreensfct.hpp" />
     <ClInclude Include="ql\experimental\finitedifferences\fdmlocalvolfwdop.hpp" />
+    <ClInclude Include="ql\experimental\finitedifferences\fdornsteinuhlenbeckvanillaengine.hpp" />
     <ClInclude Include="ql\experimental\finitedifferences\hestonrndcalculator.hpp" />
     <ClInclude Include="ql\experimental\finitedifferences\localvolrndcalculator.hpp" />
     <ClInclude Include="ql\experimental\finitedifferences\riskneutraldensitycalculator.hpp" />
@@ -497,12 +498,14 @@
     <ClInclude Include="ql\experimental\models\all.hpp" />
     <ClInclude Include="ql\experimental\models\hestonslvfdmmodel.hpp" />
     <ClInclude Include="ql\experimental\models\hestonslvmcmodel.hpp" />
+    <ClInclude Include="ql\experimental\models\normalclvmodel.hpp" />
     <ClInclude Include="ql\indexes\ibor\aonia.hpp" />
     <ClInclude Include="ql\indexes\ibor\bbsw.hpp" />
     <ClInclude Include="ql\indexes\ibor\bkbm.hpp" />
     <ClInclude Include="ql\indexes\ibor\nzocr.hpp" />
     <ClInclude Include="ql\math\polynomialmathfunction.hpp" />
     <ClInclude Include="ql\math\pascaltriangle.hpp" />
+    <ClInclude Include="ql\methods\finitedifferences\operators\fdmornsteinuhlenbeckop.hpp" />
     <ClInclude Include="ql\rebatedexercise.hpp" />
     <ClInclude Include="ql\experimental\finitedifferences\dynprogvppintrinsicvalueengine.hpp" />
     <ClInclude Include="ql\experimental\finitedifferences\fdextoujumpvanillaengine.hpp" />
@@ -843,6 +846,7 @@
     <ClInclude Include="ql\math\interpolations\interpolation2d.hpp" />
     <ClInclude Include="ql\math\interpolations\kernelinterpolation.hpp" />
     <ClInclude Include="ql\math\interpolations\kernelinterpolation2d.hpp" />
+    <ClInclude Include="ql\math\interpolations\lagrangeinterpolation.hpp" />
     <ClInclude Include="ql\math\interpolations\linearinterpolation.hpp" />
     <ClInclude Include="ql\math\interpolations\loginterpolation.hpp" />
     <ClInclude Include="ql\math\interpolations\mixedinterpolation.hpp" />
@@ -910,6 +914,7 @@
     <ClInclude Include="ql\math\randomnumbers\rngtraits.hpp" />
     <ClInclude Include="ql\math\randomnumbers\seedgenerator.hpp" />
     <ClInclude Include="ql\math\randomnumbers\sobolrsg.hpp" />
+    <ClInclude Include="ql\math\randomnumbers\stochasticcollocationinvcdf.hpp" />
     <ClInclude Include="ql\math\solvers1d\all.hpp" />
     <ClInclude Include="ql\math\solvers1d\bisection.hpp" />
     <ClInclude Include="ql\math\solvers1d\brent.hpp" />
@@ -1770,14 +1775,17 @@
     <ClCompile Include="ql\experimental\finitedifferences\bsmrndcalculator.cpp" />
     <ClCompile Include="ql\experimental\finitedifferences\fdmhestongreensfct.cpp" />
     <ClCompile Include="ql\experimental\finitedifferences\fdmlocalvolfwdop.cpp" />
+    <ClCompile Include="ql\experimental\finitedifferences\fdornsteinuhlenbeckvanillaengine.cpp" />
     <ClCompile Include="ql\experimental\finitedifferences\hestonrndcalculator.cpp" />
     <ClCompile Include="ql\experimental\finitedifferences\localvolrndcalculator.cpp" />
     <ClCompile Include="ql\experimental\finitedifferences\riskneutraldensitycalculator.cpp" />
     <ClCompile Include="ql\experimental\finitedifferences\squarerootprocessrndcalculator.cpp" />
     <ClCompile Include="ql\experimental\models\hestonslvfdmmodel.cpp" />
     <ClCompile Include="ql\experimental\models\hestonslvmcmodel.cpp" />
+    <ClCompile Include="ql\experimental\models\normalclvmodel.cpp" />
     <ClCompile Include="ql\math\polynomialmathfunction.cpp" />
     <ClCompile Include="ql\math\pascaltriangle.cpp" />
+    <ClCompile Include="ql\methods\finitedifferences\operators\fdmornsteinuhlenbeckop.cpp" />
     <ClCompile Include="ql\patterns\observable.cpp" />
     <ClCompile Include="ql\rebatedexercise.cpp" />
     <ClInclude Include="ql\experimental\finitedifferences\all.hpp" />
@@ -2031,6 +2039,7 @@
     <ClCompile Include="ql\math\randomnumbers\primitivepolynomials.cpp" />
     <ClCompile Include="ql\math\randomnumbers\seedgenerator.cpp" />
     <ClCompile Include="ql\math\randomnumbers\sobolrsg.cpp" />
+    <ClCompile Include="ql\math\randomnumbers\stochasticcollocationinvcdf.cpp" />
     <ClCompile Include="ql\math\optimization\armijo.cpp" />
     <ClCompile Include="ql\math\optimization\bfgs.cpp" />
     <ClCompile Include="ql\math\optimization\conjugategradient.cpp" />
diff --git a/QuantLib.vcxproj.filters b/QuantLib.vcxproj.filters
index 6022d9e5f30..23185e11b34 100644
--- a/QuantLib.vcxproj.filters
+++ b/QuantLib.vcxproj.filters
@@ -984,6 +984,9 @@
     <ClInclude Include="ql\math\interpolations\linearinterpolation.hpp">
       <Filter>math\interpolations</Filter>
     </ClInclude>
+    <ClInclude Include="ql\math\interpolations\lagrangeinterpolation.hpp">
+      <Filter>math\interpolations</Filter>
+    </ClInclude>
     <ClInclude Include="ql\math\interpolations\loginterpolation.hpp">
       <Filter>math\interpolations</Filter>
     </ClInclude>
@@ -1182,6 +1185,9 @@
     <ClInclude Include="ql\math\randomnumbers\sobolrsg.hpp">
       <Filter>math\randomnumbers</Filter>
     </ClInclude>
+    <ClInclude Include="ql\math\randomnumbers\stochasticcollocationinvcdf.hpp">
+      <Filter>math\randomnumbers</Filter>
+    </ClInclude>
     <ClInclude Include="ql\math\solvers1d\all.hpp">
       <Filter>math\solvers1D</Filter>
     </ClInclude>
@@ -4039,6 +4045,9 @@
     <ClInclude Include="ql\experimental\models\hestonslvmcmodel.hpp">
       <Filter>experimental\models</Filter>
     </ClInclude>
+    <ClInclude Include="ql\experimental\models\normalclvmodel.hpp">
+      <Filter>experimental\models</Filter>
+    </ClInclude>
     <ClInclude Include="ql\termstructures\volatility\equityfx\fixedlocalvolsurface.hpp">
       <Filter>termstructures\volatility\equityfx</Filter>
     </ClInclude>
@@ -4048,6 +4057,9 @@
     <ClInclude Include="ql\experimental\finitedifferences\fdmlocalvolfwdop.hpp">
       <Filter>experimental\finitedifferences</Filter>
     </ClInclude>
+    <ClInclude Include="ql\experimental\finitedifferences\fdornsteinuhlenbeckvanillaengine.hpp">
+      <Filter>experimental\finitedifferences</Filter>
+    </ClInclude>
     <ClInclude Include="ql\termstructures\volatility\equityfx\noexceptlocalvolsurface.hpp">
       <Filter>termstructures\volatility\equityfx</Filter>
     </ClInclude>
@@ -4090,6 +4102,9 @@
     <ClInclude Include="ql\termstructures\volatility\equityfx\hestonblackvolsurface.hpp">
       <Filter>termstructures\volatility\equityfx</Filter>
     </ClInclude>
+    <ClInclude Include="ql\methods\finitedifferences\operators\fdmornsteinuhlenbeckop.hpp">
+      <Filter>methods\finitedifferences\operators</Filter>
+    </ClInclude>
     <ClInclude Include="ql\indexes\ibor\aonia.hpp">
       <Filter>indexes\ibor</Filter>
     </ClInclude>
@@ -4569,6 +4584,9 @@
     <ClCompile Include="ql\math\randomnumbers\sobolrsg.cpp">
       <Filter>math\randomnumbers</Filter>
     </ClCompile>
+    <ClCompile Include="ql\math\randomnumbers\stochasticcollocationinvcdf.cpp">
+      <Filter>math\randomnumbers</Filter>
+    </ClCompile>
     <ClCompile Include="ql\math\optimization\armijo.cpp">
       <Filter>math\optimization</Filter>
     </ClCompile>
@@ -6566,6 +6584,9 @@
     <ClCompile Include="ql\experimental\finitedifferences\squarerootprocessrndcalculator.cpp">
       <Filter>experimental\finitedifferences</Filter>
     </ClCompile>
+    <ClCompile Include="ql\experimental\finitedifferences\fdornsteinuhlenbeckvanillaengine.cpp">
+      <Filter>experimental\finitedifferences</Filter>
+    </ClCompile>
     <ClCompile Include="ql\experimental\models\hestonslvfdmmodel.cpp">
       <Filter>experimental\models</Filter>
     </ClCompile>
@@ -6611,5 +6632,11 @@
     <ClCompile Include="ql\experimental\math\particleswarmoptimization.cpp">
       <Filter>experimental\math</Filter>
     </ClCompile>
+    <ClCompile Include="ql\experimental\models\normalclvmodel.cpp">
+      <Filter>experimental\models</Filter>
+    </ClCompile>
+    <ClCompile Include="ql\methods\finitedifferences\operators\fdmornsteinuhlenbeckop.cpp">
+      <Filter>methods\finitedifferences\operators</Filter>
+    </ClCompile>
   </ItemGroup>
-</Project>
\ No newline at end of file
+</Project>
diff --git a/QuantLib_vc14.sln b/QuantLib_vc14.sln
index ac21bfbe694..5bcfcc4d3eb 100644
--- a/QuantLib_vc14.sln
+++ b/QuantLib_vc14.sln
@@ -1,6 +1,6 @@
 Microsoft Visual Studio Solution File, Format Version 12.00
 # Visual Studio 14
-VisualStudioVersion = 14.0.22823.1
+VisualStudioVersion = 14.0.23107.0
 MinimumVisualStudioVersion = 10.0.40219.1
 Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "QuantLib", "QuantLib.vcxproj", "{AD0A27DA-91DA-46A2-ACBD-296C419ED3AA}"
 EndProject
@@ -83,7 +83,6 @@ Global
 		{A613045C-34AF-4706-AA3C-730C92524F74}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{A613045C-34AF-4706-AA3C-730C92524F74}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{A613045C-34AF-4706-AA3C-730C92524F74}.Release|Win32.ActiveCfg = Release|Win32
-		{A613045C-34AF-4706-AA3C-730C92524F74}.Release|Win32.Build.0 = Release|Win32
 		{A613045C-34AF-4706-AA3C-730C92524F74}.Release|x64.ActiveCfg = Release|x64
 		{A613045C-34AF-4706-AA3C-730C92524F74}.Release|x64.Build.0 = Release|x64
 		{4EAC6A0E-20F2-4B5A-8250-7E930CCE3AD0}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -99,7 +98,6 @@ Global
 		{4EAC6A0E-20F2-4B5A-8250-7E930CCE3AD0}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{4EAC6A0E-20F2-4B5A-8250-7E930CCE3AD0}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{4EAC6A0E-20F2-4B5A-8250-7E930CCE3AD0}.Release|Win32.ActiveCfg = Release|Win32
-		{4EAC6A0E-20F2-4B5A-8250-7E930CCE3AD0}.Release|Win32.Build.0 = Release|Win32
 		{4EAC6A0E-20F2-4B5A-8250-7E930CCE3AD0}.Release|x64.ActiveCfg = Release|x64
 		{4EAC6A0E-20F2-4B5A-8250-7E930CCE3AD0}.Release|x64.Build.0 = Release|x64
 		{EF6D982C-CF99-4442-B297-776DBECFAFC9}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -115,7 +113,6 @@ Global
 		{EF6D982C-CF99-4442-B297-776DBECFAFC9}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{EF6D982C-CF99-4442-B297-776DBECFAFC9}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{EF6D982C-CF99-4442-B297-776DBECFAFC9}.Release|Win32.ActiveCfg = Release|Win32
-		{EF6D982C-CF99-4442-B297-776DBECFAFC9}.Release|Win32.Build.0 = Release|Win32
 		{EF6D982C-CF99-4442-B297-776DBECFAFC9}.Release|x64.ActiveCfg = Release|x64
 		{EF6D982C-CF99-4442-B297-776DBECFAFC9}.Release|x64.Build.0 = Release|x64
 		{7C32702E-ED12-49F1-B476-656DD1EBCE66}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -131,7 +128,6 @@ Global
 		{7C32702E-ED12-49F1-B476-656DD1EBCE66}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{7C32702E-ED12-49F1-B476-656DD1EBCE66}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{7C32702E-ED12-49F1-B476-656DD1EBCE66}.Release|Win32.ActiveCfg = Release|Win32
-		{7C32702E-ED12-49F1-B476-656DD1EBCE66}.Release|Win32.Build.0 = Release|Win32
 		{7C32702E-ED12-49F1-B476-656DD1EBCE66}.Release|x64.ActiveCfg = Release|x64
 		{7C32702E-ED12-49F1-B476-656DD1EBCE66}.Release|x64.Build.0 = Release|x64
 		{B96E9E0A-99DA-4E9F-B8D0-941F46CDF634}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -147,7 +143,6 @@ Global
 		{B96E9E0A-99DA-4E9F-B8D0-941F46CDF634}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{B96E9E0A-99DA-4E9F-B8D0-941F46CDF634}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{B96E9E0A-99DA-4E9F-B8D0-941F46CDF634}.Release|Win32.ActiveCfg = Release|Win32
-		{B96E9E0A-99DA-4E9F-B8D0-941F46CDF634}.Release|Win32.Build.0 = Release|Win32
 		{B96E9E0A-99DA-4E9F-B8D0-941F46CDF634}.Release|x64.ActiveCfg = Release|x64
 		{B96E9E0A-99DA-4E9F-B8D0-941F46CDF634}.Release|x64.Build.0 = Release|x64
 		{1B660588-A923-4D84-9092-16DA67869773}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -163,7 +158,6 @@ Global
 		{1B660588-A923-4D84-9092-16DA67869773}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{1B660588-A923-4D84-9092-16DA67869773}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{1B660588-A923-4D84-9092-16DA67869773}.Release|Win32.ActiveCfg = Release|Win32
-		{1B660588-A923-4D84-9092-16DA67869773}.Release|Win32.Build.0 = Release|Win32
 		{1B660588-A923-4D84-9092-16DA67869773}.Release|x64.ActiveCfg = Release|x64
 		{1B660588-A923-4D84-9092-16DA67869773}.Release|x64.Build.0 = Release|x64
 		{7FF22935-8C7D-4903-908C-B77A9CDBA840}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -179,7 +173,6 @@ Global
 		{7FF22935-8C7D-4903-908C-B77A9CDBA840}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{7FF22935-8C7D-4903-908C-B77A9CDBA840}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{7FF22935-8C7D-4903-908C-B77A9CDBA840}.Release|Win32.ActiveCfg = Release|Win32
-		{7FF22935-8C7D-4903-908C-B77A9CDBA840}.Release|Win32.Build.0 = Release|Win32
 		{7FF22935-8C7D-4903-908C-B77A9CDBA840}.Release|x64.ActiveCfg = Release|x64
 		{7FF22935-8C7D-4903-908C-B77A9CDBA840}.Release|x64.Build.0 = Release|x64
 		{940A0AFC-9F9F-4797-A0FF-99543F67C1D9}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -195,7 +188,6 @@ Global
 		{940A0AFC-9F9F-4797-A0FF-99543F67C1D9}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{940A0AFC-9F9F-4797-A0FF-99543F67C1D9}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{940A0AFC-9F9F-4797-A0FF-99543F67C1D9}.Release|Win32.ActiveCfg = Release|Win32
-		{940A0AFC-9F9F-4797-A0FF-99543F67C1D9}.Release|Win32.Build.0 = Release|Win32
 		{940A0AFC-9F9F-4797-A0FF-99543F67C1D9}.Release|x64.ActiveCfg = Release|x64
 		{940A0AFC-9F9F-4797-A0FF-99543F67C1D9}.Release|x64.Build.0 = Release|x64
 		{C3E22CAD-0CAF-42DC-ADE0-B2FF4F644BCE}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -211,7 +203,6 @@ Global
 		{C3E22CAD-0CAF-42DC-ADE0-B2FF4F644BCE}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{C3E22CAD-0CAF-42DC-ADE0-B2FF4F644BCE}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{C3E22CAD-0CAF-42DC-ADE0-B2FF4F644BCE}.Release|Win32.ActiveCfg = Release|Win32
-		{C3E22CAD-0CAF-42DC-ADE0-B2FF4F644BCE}.Release|Win32.Build.0 = Release|Win32
 		{C3E22CAD-0CAF-42DC-ADE0-B2FF4F644BCE}.Release|x64.ActiveCfg = Release|x64
 		{C3E22CAD-0CAF-42DC-ADE0-B2FF4F644BCE}.Release|x64.Build.0 = Release|x64
 		{4E262A25-90B4-449A-BFC0-95311CADF91D}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -227,7 +218,6 @@ Global
 		{4E262A25-90B4-449A-BFC0-95311CADF91D}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{4E262A25-90B4-449A-BFC0-95311CADF91D}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{4E262A25-90B4-449A-BFC0-95311CADF91D}.Release|Win32.ActiveCfg = Release|Win32
-		{4E262A25-90B4-449A-BFC0-95311CADF91D}.Release|Win32.Build.0 = Release|Win32
 		{4E262A25-90B4-449A-BFC0-95311CADF91D}.Release|x64.ActiveCfg = Release|x64
 		{4E262A25-90B4-449A-BFC0-95311CADF91D}.Release|x64.Build.0 = Release|x64
 		{0214688B-CE8A-4446-9BDD-1AE7F486EB8B}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -243,7 +233,6 @@ Global
 		{0214688B-CE8A-4446-9BDD-1AE7F486EB8B}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{0214688B-CE8A-4446-9BDD-1AE7F486EB8B}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{0214688B-CE8A-4446-9BDD-1AE7F486EB8B}.Release|Win32.ActiveCfg = Release|Win32
-		{0214688B-CE8A-4446-9BDD-1AE7F486EB8B}.Release|Win32.Build.0 = Release|Win32
 		{0214688B-CE8A-4446-9BDD-1AE7F486EB8B}.Release|x64.ActiveCfg = Release|x64
 		{0214688B-CE8A-4446-9BDD-1AE7F486EB8B}.Release|x64.Build.0 = Release|x64
 		{65F5530B-D97E-4BDB-949F-9C31C56104B0}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -259,7 +248,6 @@ Global
 		{65F5530B-D97E-4BDB-949F-9C31C56104B0}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{65F5530B-D97E-4BDB-949F-9C31C56104B0}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{65F5530B-D97E-4BDB-949F-9C31C56104B0}.Release|Win32.ActiveCfg = Release|Win32
-		{65F5530B-D97E-4BDB-949F-9C31C56104B0}.Release|Win32.Build.0 = Release|Win32
 		{65F5530B-D97E-4BDB-949F-9C31C56104B0}.Release|x64.ActiveCfg = Release|x64
 		{65F5530B-D97E-4BDB-949F-9C31C56104B0}.Release|x64.Build.0 = Release|x64
 		{C8AFC9D1-704C-4AB5-911B-3C93C27C63D0}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -275,7 +263,6 @@ Global
 		{C8AFC9D1-704C-4AB5-911B-3C93C27C63D0}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{C8AFC9D1-704C-4AB5-911B-3C93C27C63D0}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{C8AFC9D1-704C-4AB5-911B-3C93C27C63D0}.Release|Win32.ActiveCfg = Release|Win32
-		{C8AFC9D1-704C-4AB5-911B-3C93C27C63D0}.Release|Win32.Build.0 = Release|Win32
 		{C8AFC9D1-704C-4AB5-911B-3C93C27C63D0}.Release|x64.ActiveCfg = Release|x64
 		{C8AFC9D1-704C-4AB5-911B-3C93C27C63D0}.Release|x64.Build.0 = Release|x64
 		{47CE2A41-091A-42A3-B40D-F6F0DD689349}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -291,7 +278,6 @@ Global
 		{47CE2A41-091A-42A3-B40D-F6F0DD689349}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{47CE2A41-091A-42A3-B40D-F6F0DD689349}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{47CE2A41-091A-42A3-B40D-F6F0DD689349}.Release|Win32.ActiveCfg = Release|Win32
-		{47CE2A41-091A-42A3-B40D-F6F0DD689349}.Release|Win32.Build.0 = Release|Win32
 		{47CE2A41-091A-42A3-B40D-F6F0DD689349}.Release|x64.ActiveCfg = Release|x64
 		{47CE2A41-091A-42A3-B40D-F6F0DD689349}.Release|x64.Build.0 = Release|x64
 		{B62FC7BE-C1BC-4AB0-BB11-BBF627CC3BA6}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -307,7 +293,6 @@ Global
 		{B62FC7BE-C1BC-4AB0-BB11-BBF627CC3BA6}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{B62FC7BE-C1BC-4AB0-BB11-BBF627CC3BA6}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{B62FC7BE-C1BC-4AB0-BB11-BBF627CC3BA6}.Release|Win32.ActiveCfg = Release|Win32
-		{B62FC7BE-C1BC-4AB0-BB11-BBF627CC3BA6}.Release|Win32.Build.0 = Release|Win32
 		{B62FC7BE-C1BC-4AB0-BB11-BBF627CC3BA6}.Release|x64.ActiveCfg = Release|x64
 		{B62FC7BE-C1BC-4AB0-BB11-BBF627CC3BA6}.Release|x64.Build.0 = Release|x64
 		{D3E10D43-57DF-42D7-AEA6-44AD48CAE119}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -323,7 +308,6 @@ Global
 		{D3E10D43-57DF-42D7-AEA6-44AD48CAE119}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{D3E10D43-57DF-42D7-AEA6-44AD48CAE119}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{D3E10D43-57DF-42D7-AEA6-44AD48CAE119}.Release|Win32.ActiveCfg = Release|Win32
-		{D3E10D43-57DF-42D7-AEA6-44AD48CAE119}.Release|Win32.Build.0 = Release|Win32
 		{D3E10D43-57DF-42D7-AEA6-44AD48CAE119}.Release|x64.ActiveCfg = Release|x64
 		{D3E10D43-57DF-42D7-AEA6-44AD48CAE119}.Release|x64.Build.0 = Release|x64
 		{E4098DFF-FD34-48A6-955F-C1E3F97FBA26}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -339,7 +323,6 @@ Global
 		{E4098DFF-FD34-48A6-955F-C1E3F97FBA26}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{E4098DFF-FD34-48A6-955F-C1E3F97FBA26}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{E4098DFF-FD34-48A6-955F-C1E3F97FBA26}.Release|Win32.ActiveCfg = Release|Win32
-		{E4098DFF-FD34-48A6-955F-C1E3F97FBA26}.Release|Win32.Build.0 = Release|Win32
 		{E4098DFF-FD34-48A6-955F-C1E3F97FBA26}.Release|x64.ActiveCfg = Release|x64
 		{E4098DFF-FD34-48A6-955F-C1E3F97FBA26}.Release|x64.Build.0 = Release|x64
 		{43A17E5B-EC94-4EB5-9D68-788BF234AE1F}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -355,7 +338,6 @@ Global
 		{43A17E5B-EC94-4EB5-9D68-788BF234AE1F}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{43A17E5B-EC94-4EB5-9D68-788BF234AE1F}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{43A17E5B-EC94-4EB5-9D68-788BF234AE1F}.Release|Win32.ActiveCfg = Release|Win32
-		{43A17E5B-EC94-4EB5-9D68-788BF234AE1F}.Release|Win32.Build.0 = Release|Win32
 		{43A17E5B-EC94-4EB5-9D68-788BF234AE1F}.Release|x64.ActiveCfg = Release|x64
 		{43A17E5B-EC94-4EB5-9D68-788BF234AE1F}.Release|x64.Build.0 = Release|x64
 		{D85EA161-361A-4CAC-8B77-39E174FCFFE2}.Debug (static runtime)|Win32.ActiveCfg = Debug (static runtime)|Win32
@@ -371,7 +353,6 @@ Global
 		{D85EA161-361A-4CAC-8B77-39E174FCFFE2}.Release (static runtime)|x64.ActiveCfg = Release (static runtime)|x64
 		{D85EA161-361A-4CAC-8B77-39E174FCFFE2}.Release (static runtime)|x64.Build.0 = Release (static runtime)|x64
 		{D85EA161-361A-4CAC-8B77-39E174FCFFE2}.Release|Win32.ActiveCfg = Release|Win32
-		{D85EA161-361A-4CAC-8B77-39E174FCFFE2}.Release|Win32.Build.0 = Release|Win32
 		{D85EA161-361A-4CAC-8B77-39E174FCFFE2}.Release|x64.ActiveCfg = Release|x64
 		{D85EA161-361A-4CAC-8B77-39E174FCFFE2}.Release|x64.Build.0 = Release|x64
 	EndGlobalSection
diff --git a/ql/experimental/finitedifferences/Makefile.am b/ql/experimental/finitedifferences/Makefile.am
index 1ce3f470f3b..adb45e05dff 100644
--- a/ql/experimental/finitedifferences/Makefile.am
+++ b/ql/experimental/finitedifferences/Makefile.am
@@ -29,6 +29,7 @@ this_include_HEADERS = \
 	fdmvppstepcondition.hpp \
 	fdmvppstepconditionfactory.hpp \
 	fdmzabrop.hpp \
+	fdornsteinuhlenbeckvanillaengine.hpp \
 	fdsimpleextoujumpswingengine.hpp \
     fdsimpleextoustorageengine.hpp \
 	fdsimpleklugeextouvppengine.hpp \
@@ -60,6 +61,7 @@ libMultiDimFDM_la_SOURCES = \
 	fdmvppstepcondition.cpp \
 	fdmvppstepconditionfactory.cpp \
 	fdmzabrop.cpp \
+	fdornsteinuhlenbeckvanillaengine.cpp \
 	fdsimpleextoujumpswingengine.cpp \
     fdsimpleextoustorageengine.cpp \
     fdsimpleklugeextouvppengine.cpp \
diff --git a/ql/experimental/finitedifferences/all.hpp b/ql/experimental/finitedifferences/all.hpp
index 9756dad2bec..d11c1061667 100644
--- a/ql/experimental/finitedifferences/all.hpp
+++ b/ql/experimental/finitedifferences/all.hpp
@@ -26,6 +26,7 @@
 #include <ql/experimental/finitedifferences/fdmvppstepcondition.hpp>
 #include <ql/experimental/finitedifferences/fdmvppstepconditionfactory.hpp>
 #include <ql/experimental/finitedifferences/fdmzabrop.hpp>
+#include <ql/experimental/finitedifferences/fdornsteinuhlenbeckvanillaengine.hpp>
 #include <ql/experimental/finitedifferences/fdsimpleextoujumpswingengine.hpp>
 #include <ql/experimental/finitedifferences/fdsimpleextoustorageengine.hpp>
 #include <ql/experimental/finitedifferences/fdsimpleklugeextouvppengine.hpp>
diff --git a/ql/experimental/finitedifferences/fdornsteinuhlenbeckvanillaengine.cpp b/ql/experimental/finitedifferences/fdornsteinuhlenbeckvanillaengine.cpp
new file mode 100644
index 00000000000..cda4590b4b6
--- /dev/null
+++ b/ql/experimental/finitedifferences/fdornsteinuhlenbeckvanillaengine.cpp
@@ -0,0 +1,125 @@
+/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+
+/*
+ Copyright (C) 2016 Klaus Spanderen
+
+ This file is part of QuantLib, a free-software/open-source library
+ for financial quantitative analysts and developers - http://quantlib.org/
+
+ QuantLib is free software: you can redistribute it and/or modify it
+ under the terms of the QuantLib license.  You should have received a
+ copy of the license along with this program; if not, please email
+ <quantlib-dev@lists.sf.net>. The license is also available online at
+ <http://quantlib.org/license.shtml>.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE.  See the license for more details.
+*/
+
+#include <ql/exercise.hpp>
+#include <ql/termstructures/yieldtermstructure.hpp>
+#include <ql/processes/ornsteinuhlenbeckprocess.hpp>
+#include <ql/methods/finitedifferences/solvers/fdm1dimsolver.hpp>
+#include <ql/methods/finitedifferences/meshers/fdmmeshercomposite.hpp>
+#include <ql/methods/finitedifferences/meshers/fdmsimpleprocess1dmesher.hpp>
+#include <ql/methods/finitedifferences/operators/fdmlinearoplayout.hpp>
+#include <ql/methods/finitedifferences/operators/fdmornsteinuhlenbeckop.hpp>
+#include <ql/methods/finitedifferences/utilities/fdminnervaluecalculator.hpp>
+#include <ql/methods/finitedifferences/stepconditions/fdmstepconditioncomposite.hpp>
+#include <ql/experimental/finitedifferences/fdornsteinuhlenbeckvanillaengine.hpp>
+
+namespace QuantLib {
+
+    namespace {
+        class FdmOUInnerValue : public FdmInnerValueCalculator {
+          public:
+            FdmOUInnerValue(
+                const boost::shared_ptr<Payoff>& payoff,
+                const boost::shared_ptr<FdmMesher>& mesher,
+                Size direction)
+          : payoff_(payoff), mesher_(mesher), direction_ (direction) { }
+
+
+            Real innerValue(const FdmLinearOpIterator& iter, Time t) {
+                const Real s = mesher_->location(iter, direction_);
+                return payoff_->operator()(s);
+            }
+
+            Real avgInnerValue(const FdmLinearOpIterator& iter, Time t) {
+                return innerValue(iter, t);
+            }
+
+          private:
+            const boost::shared_ptr<Payoff> payoff_;
+            const boost::shared_ptr<FdmMesher> mesher_;
+            const Size direction_;
+        };
+    }
+
+    FdOrnsteinUhlenbeckVanillaEngine::FdOrnsteinUhlenbeckVanillaEngine(
+            const boost::shared_ptr<OrnsteinUhlenbeckProcess>& process,
+            const boost::shared_ptr<YieldTermStructure>& rTS,
+            Size tGrid, Size xGrid, Size dampingSteps,
+            Real epsilon,
+            const FdmSchemeDesc& schemeDesc)
+    : process_(process),
+      rTS_(rTS),
+      tGrid_(tGrid), xGrid_(xGrid), dampingSteps_(dampingSteps),
+      epsilon_(epsilon),
+      schemeDesc_(schemeDesc) {
+        registerWith(process_);
+        registerWith(rTS);
+    }
+
+    void FdOrnsteinUhlenbeckVanillaEngine::calculate() const {
+
+        // 1. Mesher
+        const boost::shared_ptr<StrikedTypePayoff> payoff =
+            boost::dynamic_pointer_cast<StrikedTypePayoff>(arguments_.payoff);
+
+        const DayCounter dc = rTS_->dayCounter();
+        const Date referenceDate = rTS_->referenceDate();
+
+        const Time maturity = dc.yearFraction(
+            referenceDate, arguments_.exercise->lastDate());
+
+        const boost::shared_ptr<Fdm1dMesher> equityMesher(
+            new FdmSimpleProcess1dMesher(
+                xGrid_, process_, maturity, 1, epsilon_));
+
+        const boost::shared_ptr<FdmMesher> mesher (
+            new FdmMesherComposite(equityMesher));
+
+        // 2. Calculator
+        const boost::shared_ptr<FdmInnerValueCalculator> calculator(
+            new FdmOUInnerValue(payoff, mesher, 0));
+
+        // 3. Step conditions
+        const boost::shared_ptr<FdmStepConditionComposite> conditions =
+            FdmStepConditionComposite::vanillaComposite(
+                                    arguments_.cashFlow, arguments_.exercise,
+                                    mesher, calculator,
+                                    referenceDate, dc);
+
+        // 4. Boundary conditions
+        const FdmBoundaryConditionSet boundaries;
+
+        // 5. Solver
+        FdmSolverDesc solverDesc = { mesher, boundaries, conditions, calculator,
+                                     maturity, tGrid_, dampingSteps_ };
+
+        const boost::shared_ptr<FdmOrnsteinUhlenbackOp> op(
+            new FdmOrnsteinUhlenbackOp(mesher, process_, rTS_, boundaries, 0));
+
+        const boost::shared_ptr<Fdm1DimSolver> solver(
+                new Fdm1DimSolver(solverDesc, schemeDesc_, op));
+
+        const Real spot = process_->x0();
+
+        results_.value = solver->interpolateAt(spot);
+        results_.delta = solver->derivativeX(spot);
+        results_.gamma = solver->derivativeXX(spot);
+        results_.theta = solver->thetaAt(spot);
+    }
+}
diff --git a/ql/experimental/finitedifferences/fdornsteinuhlenbeckvanillaengine.hpp b/ql/experimental/finitedifferences/fdornsteinuhlenbeckvanillaengine.hpp
new file mode 100644
index 00000000000..dc57f05f9d0
--- /dev/null
+++ b/ql/experimental/finitedifferences/fdornsteinuhlenbeckvanillaengine.hpp
@@ -0,0 +1,58 @@
+/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+
+/*
+ Copyright (C) 2016 Klaus Spanderen
+
+ This file is part of QuantLib, a free-software/open-source library
+ for financial quantitative analysts and developers - http://quantlib.org/
+
+ QuantLib is free software: you can redistribute it and/or modify it
+ under the terms of the QuantLib license.  You should have received a
+ copy of the license along with this program; if not, please email
+ <quantlib-dev@lists.sf.net>. The license is also available online at
+ <http://quantlib.org/license.shtml>.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE.  See the license for more details.
+*/
+
+/*! \file fdblackscholesvanillaengine.hpp
+    \brief Finite-Differences Black Scholes vanilla option engine
+*/
+
+#ifndef quantlib_fd_ornstein_uhlenbeck_vanilla_engine_hpp
+#define quantlib_fd_ornstein_uhlenbeck_vanilla_engine_hpp
+
+#include <ql/pricingengine.hpp>
+#include <ql/instruments/dividendvanillaoption.hpp>
+#include <ql/methods/finitedifferences/solvers/fdmbackwardsolver.hpp>
+
+namespace QuantLib {
+
+    class YieldTermStructure;
+    class OrnsteinUhlenbeckProcess;
+
+    class FdOrnsteinUhlenbeckVanillaEngine
+         : public DividendVanillaOption::engine {
+      public:
+        // Constructor
+        FdOrnsteinUhlenbeckVanillaEngine(
+            const boost::shared_ptr<OrnsteinUhlenbeckProcess>&,
+            const boost::shared_ptr<YieldTermStructure>& rTS,
+            Size tGrid = 100, Size xGrid = 100, Size dampingSteps = 0,
+            Real epsilon = 0.0001,
+            const FdmSchemeDesc& schemeDesc = FdmSchemeDesc::Douglas());
+
+        void calculate() const;
+
+      private:
+        const boost::shared_ptr<OrnsteinUhlenbeckProcess> process_;
+        const boost::shared_ptr<YieldTermStructure> rTS_;
+        const Size tGrid_, xGrid_, dampingSteps_;
+        const Real epsilon_;
+        const FdmSchemeDesc schemeDesc_;
+    };
+}
+
+#endif
diff --git a/ql/experimental/models/Makefile.am b/ql/experimental/models/Makefile.am
index 908bf285243..9f82e08975a 100644
--- a/ql/experimental/models/Makefile.am
+++ b/ql/experimental/models/Makefile.am
@@ -5,11 +5,13 @@ this_includedir=${includedir}/${subdir}
 this_include_HEADERS = \
     all.hpp \
     hestonslvfdmmodel.hpp \
-    hestonslvmcmodel.hpp
+    hestonslvmcmodel.hpp \
+    normalclvmodel.hpp
 
 libModels_la_SOURCES = \
     hestonslvfdmmodel.cpp \
-    hestonslvmcmodel.cpp
+    hestonslvmcmodel.cpp \
+    normalclvmodel.cpp
 
 noinst_LTLIBRARIES = libModels.la
 
diff --git a/ql/experimental/models/all.hpp b/ql/experimental/models/all.hpp
index 018f8a7fe42..e72242b050c 100644
--- a/ql/experimental/models/all.hpp
+++ b/ql/experimental/models/all.hpp
@@ -3,4 +3,5 @@
 
 #include <ql/experimental/models/hestonslvfdmmodel.hpp>
 #include <ql/experimental/models/hestonslvmcmodel.hpp>
+#include <ql/experimental/models/normalclvmodel.hpp>
 
diff --git a/ql/experimental/models/hestonslvfdmmodel.cpp b/ql/experimental/models/hestonslvfdmmodel.cpp
index d75e14275b8..684e8a50411 100644
--- a/ql/experimental/models/hestonslvfdmmodel.cpp
+++ b/ql/experimental/models/hestonslvfdmmodel.cpp
@@ -281,10 +281,6 @@ namespace QuantLib {
         registerWith(hestonModel_);
     }
 
-    void HestonSLVFDMModel::update() {
-        notifyObservers();
-    }
-
     boost::shared_ptr<HestonProcess> HestonSLVFDMModel::hestonProcess() const {
         return hestonModel_->process();
     }
diff --git a/ql/experimental/models/hestonslvfdmmodel.hpp b/ql/experimental/models/hestonslvfdmmodel.hpp
index 5ec07a560ca..06f1ee714d4 100644
--- a/ql/experimental/models/hestonslvfdmmodel.hpp
+++ b/ql/experimental/models/hestonslvfdmmodel.hpp
@@ -77,8 +77,6 @@ class SimpleQuote;
             const bool logging = false,
             const std::vector<Date>& mandatoryDates = std::vector<Date>());
 
-        void update();
-
         boost::shared_ptr<HestonProcess> hestonProcess() const;
         boost::shared_ptr<LocalVolTermStructure> localVol() const;
         boost::shared_ptr<LocalVolTermStructure> leverageFunction() const;
diff --git a/ql/experimental/models/hestonslvmcmodel.cpp b/ql/experimental/models/hestonslvmcmodel.cpp
index 5f695592121..e01e1057f30 100644
--- a/ql/experimental/models/hestonslvmcmodel.cpp
+++ b/ql/experimental/models/hestonslvmcmodel.cpp
@@ -71,10 +71,6 @@ namespace QuantLib {
                 std::max(Size(2), Size(gridTimes.back()*timeStepsPerYear)));
     }
 
-    void HestonSLVMCModel::update() {
-        notifyObservers();
-    }
-
     boost::shared_ptr<HestonProcess> HestonSLVMCModel::hestonProcess() const {
         return hestonModel_->process();
     }
diff --git a/ql/experimental/models/hestonslvmcmodel.hpp b/ql/experimental/models/hestonslvmcmodel.hpp
index 8e87d75100b..5203e15510f 100644
--- a/ql/experimental/models/hestonslvmcmodel.hpp
+++ b/ql/experimental/models/hestonslvmcmodel.hpp
@@ -52,8 +52,6 @@ namespace QuantLib {
             Size calibrationPaths = (1 << 15),
             const std::vector<Date>& mandatoryDates = std::vector<Date>());
 
-        void update();
-
         boost::shared_ptr<HestonProcess> hestonProcess() const;
         boost::shared_ptr<LocalVolTermStructure> localVol() const;
         boost::shared_ptr<LocalVolTermStructure> leverageFunction() const;
diff --git a/ql/experimental/models/normalclvmodel.cpp b/ql/experimental/models/normalclvmodel.cpp
new file mode 100644
index 00000000000..92065fd95a0
--- /dev/null
+++ b/ql/experimental/models/normalclvmodel.cpp
@@ -0,0 +1,201 @@
+/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+
+/*
+ Copyright (C) 2016 Klaus Spanderen
+
+ This file is part of QuantLib, a free-software/open-source library
+ for financial quantitative analysts and developers - http://quantlib.org/
+
+ QuantLib is free software: you can redistribute it and/or modify it
+ under the terms of the QuantLib license.  You should have received a
+ copy of the license along with this program; if not, please email
+ <quantlib-dev@lists.sf.net>. The license is also available online at
+ <http://quantlib.org/license.shtml>.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE.  See the license for more details.
+*/
+
+/*! \file normalclvmodel.cpp
+*/
+
+#include <ql/instruments/vanillaoption.hpp>
+#include <ql/math/functional.hpp>
+#include <ql/math/solvers1d/brent.hpp>
+#include <ql/math/integrals/gaussianquadratures.hpp>
+#include <ql/math/distributions/normaldistribution.hpp>
+#include <ql/math/interpolations/linearinterpolation.hpp>
+#include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp>
+#include <ql/experimental/models/normalclvmodel.hpp>
+
+#include <boost/bind.hpp>
+#include <boost/make_shared.hpp>
+
+namespace QuantLib {
+
+    NormalCLVModel::NormalCLVModel(
+        const boost::shared_ptr<GeneralizedBlackScholesProcess>& bsProcess,
+        const boost::shared_ptr<OrnsteinUhlenbeckProcess>& ouProcess,
+        const std::vector<Date>& maturityDates,
+        Size lagrangeOrder, Real pMax, Real pMin)
+    : x_(M_SQRT2*GaussHermiteIntegration(lagrangeOrder).x()),
+      sigma_( (pMax != Null<Real>())
+              ? x_.back() / InverseCumulativeNormal()(pMax)
+              : (pMin != Null<Real>())
+                  ? x_.front() / InverseCumulativeNormal()(pMin)
+                  : 1.0),
+      bsProcess_    (bsProcess),
+      ouProcess_    (ouProcess),
+      maturityDates_(maturityDates),
+      pricingEngine_(boost::make_shared<AnalyticEuropeanEngine>(bsProcess_)) {
+
+        registerWith(bsProcess_);
+        registerWith(ouProcess_);
+
+        maturityTimes_.reserve(maturityDates.size());
+        std::transform(maturityDates_.begin(), maturityDates_.end(),
+            std::back_inserter(maturityTimes_),
+            boost::bind(&GeneralizedBlackScholesProcess::time,
+                        bsProcess_, _1));
+        std::sort(maturityTimes_.begin(), maturityTimes_.end());
+    }
+
+    Real NormalCLVModel::cdf(const Date& d, Real k) const {
+        //k = std::max(k, 1e-3*bsProcess_->x0());
+
+        const Handle<BlackVolTermStructure> volTS
+            = bsProcess_->blackVolatility();
+
+        const Real dk = 1e-3*k;
+        const Real dvol_dk
+            = (volTS->blackVol(d, k+dk) - volTS->blackVol(d, k-dk)) / (2*dk);
+
+        const DiscountFactor df
+            = bsProcess_->riskFreeRate()->discount(d);
+
+        const boost::shared_ptr<Exercise> exercise
+            = boost::make_shared<EuropeanExercise>(d);
+
+        if (bsProcess_->x0() <=k) {
+            VanillaOption option(
+                boost::make_shared<PlainVanillaPayoff>(Option::Call, k),
+                exercise);
+            option.setPricingEngine(pricingEngine_);
+
+            return 1.0 + (  option.strikeSensitivity()
+                          + option.vega() * dvol_dk) /df;
+        }
+        else {
+            VanillaOption option(
+                boost::make_shared<PlainVanillaPayoff>(Option::Put, k),
+                exercise);
+            option.setPricingEngine(pricingEngine_);
+
+            return (  option.strikeSensitivity()
+                    + option.vega() * dvol_dk) /df;
+
+        }
+    }
+
+
+    Real NormalCLVModel::invCDF(const Date& d, Real q) const {
+        const Real fwd = bsProcess_->x0()
+            / bsProcess_->riskFreeRate()->discount(d, true)
+            * bsProcess_->dividendYield()->discount(d, true);
+
+        const Volatility atmVol =
+            bsProcess_->blackVolatility()->blackVol(d, fwd, true);
+
+        const Real atmX = InverseCumulativeNormal()(q);
+        const Time t = bsProcess_->time(d);
+
+        const Real guess = fwd*std::exp(atmVol*atmX*std::sqrt(t));
+
+
+        Real lower = guess;
+        while (guess/lower < 65535.0 && cdf(d, lower) > q)
+            lower*=0.5;
+
+        Real upper = guess;
+        while (upper/guess < 65535.0 && cdf(d, upper) < q) upper*=2;
+
+        QL_REQUIRE(guess/lower < 65535.0 && upper/guess < 65535.0,
+                "Could not find an start interval with ("
+                << lower << ", " << upper << ") -> ("
+                << cdf(d, lower) << ", " << cdf(d, upper) << ")");
+
+        return Brent().solve(
+            compose(std::bind2nd(std::minus<Real>(), q),
+                boost::function<Real(Real)>(
+                    boost::bind(&NormalCLVModel::cdf, this, d, _1))),
+            1e-10, 0.5*(lower+upper), lower, upper);
+    }
+
+    Disposable<Array> NormalCLVModel::collocationPointsX(const Date& d) const {
+        const Time t = bsProcess_->time(d);
+
+        const Real expectation
+            = ouProcess_->expectation(0.0, ouProcess_->x0(), t);
+        const Real stdDeviation
+            = ouProcess_->stdDeviation(0.0, ouProcess_->x0(), t);
+
+        return expectation + stdDeviation*x_;
+    }
+
+    Disposable<Array> NormalCLVModel::collocationPointsY(const Date& d) const {
+        Array s(x_.size());
+
+        CumulativeNormalDistribution N;
+        for (Size i=0, n=s.size(); i < n; ++i) {
+            s[i] = invCDF(d, N(x_[i]/sigma_));
+        }
+
+        return s;
+    }
+
+
+    boost::function<Real(Time, Real)> NormalCLVModel::g() const {
+        calculate();
+        return g_;
+    }
+
+    NormalCLVModel::MappingFunction::MappingFunction(
+        const NormalCLVModel& model)
+    : y_(model.x_.size()),
+      sigma_(model.sigma_),
+      ouProcess_(model.ouProcess_),
+      data_(boost::make_shared<InterpolationData>(model)) {
+
+        for (Size i=0; i < data_->s_.columns(); ++i) {
+            const Array y = model.collocationPointsY(model.maturityDates_[i]);
+            std::copy(y.begin(), y.end(), data_->s_.column_begin(i));
+        }
+
+        for (Size i=0; i < data_->s_.rows(); ++i) {
+            data_->interpl_.push_back(
+                LinearInterpolation(data_->t_.begin(), data_->t_.end(),
+                                    data_->s_.row_begin(i)));
+        }
+    }
+
+
+    Real NormalCLVModel::MappingFunction::operator()(Time t, Real x) const {
+        for (Size i=0; i < y_.size(); ++i) {
+            y_[i] = data_->interpl_[i](t, true);
+        }
+
+        const Real expectation
+            = ouProcess_->expectation(0.0, ouProcess_->x0(), t);
+        const Real stdDeviation
+            = ouProcess_->stdDeviation(0.0, ouProcess_->x0(), t);
+
+        const Real r = sigma_*(x-expectation)/stdDeviation;
+
+        return data_->lagrangeInterpl_.value(y_, r);
+    }
+
+    void NormalCLVModel::performCalculations() const {
+        g_ = boost::function<Real(Time, Real)>(MappingFunction(*this));
+    }
+}
diff --git a/ql/experimental/models/normalclvmodel.hpp b/ql/experimental/models/normalclvmodel.hpp
new file mode 100644
index 00000000000..c7df1c74365
--- /dev/null
+++ b/ql/experimental/models/normalclvmodel.hpp
@@ -0,0 +1,120 @@
+/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+
+/*
+ Copyright (C) 2016 Klaus Spanderen
+
+ This file is part of QuantLib, a free-software/open-source library
+ for financial quantitative analysts and developers - http://quantlib.org/
+
+ QuantLib is free software: you can redistribute it and/or modify it
+ under the terms of the QuantLib license.  You should have received a
+ copy of the license along with this program; if not, please email
+ <quantlib-dev@lists.sf.net>. The license is also available online at
+ <http://quantlib.org/license.shtml>.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE.  See the license for more details.
+*/
+
+/*! \file normalclvmodel.hpp
+    \brief CLV model with a normally distributed kernel process
+*/
+
+#ifndef quantlib_normal_clv_mc_model_hpp
+#define quantlib_normal_clv_mc_model_hpp
+
+#include <ql/handle.hpp>
+#include <ql/exercise.hpp>
+#include <ql/patterns/lazyobject.hpp>
+#include <ql/processes/blackscholesprocess.hpp>
+#include <ql/processes/ornsteinuhlenbeckprocess.hpp>
+#include <ql/math/interpolations/linearinterpolation.hpp>
+#include <ql/math/interpolations/lagrangeinterpolation.hpp>
+
+#include <boost/function.hpp>
+#include <functional>
+
+namespace QuantLib {
+    /*! References:
+
+        A. Grzelak, 2016, The CLV Framework -
+        A Fresh Look at Efficient Pricing with Smile
+
+        http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2747541
+    */
+
+    class PricingEngine;
+
+    class NormalCLVModel : public LazyObject {
+      public:
+        NormalCLVModel(
+            const boost::shared_ptr<GeneralizedBlackScholesProcess>& bsProcess,
+            const boost::shared_ptr<OrnsteinUhlenbeckProcess>& ouProcess,
+            const std::vector<Date>& maturityDates,
+            Size lagrangeOrder,
+            Real pMax = Null<Real>(),
+            Real pMin = Null<Real>());
+
+        // cumulative distribution function of the BS process
+        Real cdf(const Date& d, Real x) const;
+
+        // inverse cumulative distribution function of the BS process
+        Real invCDF(const Date& d, Real q) const;
+
+        // collocation points of the Ornstein-Uhlenbeck process
+        Disposable<Array> collocationPointsX(const Date& d) const;
+
+        // collocation points for the underlying Y
+        Disposable<Array> collocationPointsY(const Date& d) const;
+
+        // CLV mapping function
+        boost::function<Real(Time, Real)> g() const;
+
+      protected:
+        void performCalculations() const;
+
+      private:
+        class MappingFunction : public std::binary_function<Time, Real, Real> {
+          public:
+            MappingFunction(const NormalCLVModel& model);
+
+            Real operator()(Time t, Real x) const;
+
+          private:
+            mutable Array y_;
+            const Volatility sigma_;
+            const boost::shared_ptr<OrnsteinUhlenbeckProcess> ouProcess_;
+
+            struct InterpolationData {
+                InterpolationData(const NormalCLVModel& model)
+                : s_(model.x_.size(), model.maturityDates_.size()),
+                  x_(model.x_),
+                  t_(model.maturityTimes_),
+                  lagrangeInterpl_(x_.begin(), x_.end(), x_.begin()) {}
+
+                Matrix s_;
+                std::vector<LinearInterpolation> interpl_;
+
+                const Array x_;
+                const std::vector<Time> t_;
+                const LagrangeInterpolation lagrangeInterpl_;
+            };
+
+            const boost::shared_ptr<InterpolationData> data_;
+        };
+
+
+        const Array x_;
+        const Volatility sigma_;
+        const boost::shared_ptr<GeneralizedBlackScholesProcess> bsProcess_;
+        const boost::shared_ptr<OrnsteinUhlenbeckProcess> ouProcess_;
+        const std::vector<Date> maturityDates_;
+        const boost::shared_ptr<PricingEngine> pricingEngine_;
+
+        std::vector<Time> maturityTimes_;
+        mutable boost::function<Real(Time, Real)> g_;
+    };
+}
+
+#endif
diff --git a/ql/math/interpolation.hpp b/ql/math/interpolation.hpp
index b780326b6b2..d0039b63e6d 100644
--- a/ql/math/interpolation.hpp
+++ b/ql/math/interpolation.hpp
@@ -70,7 +70,7 @@ namespace QuantLib {
                 QL_REQUIRE(static_cast<int>(xEnd_-xBegin_) >= requiredPoints,
                            "not enough points to interpolate: at least " <<
                            requiredPoints <<
-                           "required, " << static_cast<int>(xEnd_-xBegin_)<< " provided");
+                           " required, " << static_cast<int>(xEnd_-xBegin_)<< " provided");
             }
             Real xMin() const {
                 return *xBegin_;
diff --git a/ql/math/interpolations/Makefile.am b/ql/math/interpolations/Makefile.am
index 98b052f096b..9a8660a7bc2 100644
--- a/ql/math/interpolations/Makefile.am
+++ b/ql/math/interpolations/Makefile.am
@@ -17,6 +17,7 @@ this_include_HEADERS = \
 	interpolation2d.hpp \
 	kernelinterpolation.hpp \
 	kernelinterpolation2d.hpp \
+	lagrangeinterpolation.hpp \
 	linearinterpolation.hpp \
 	loginterpolation.hpp \
 	mixedinterpolation.hpp \
diff --git a/ql/math/interpolations/all.hpp b/ql/math/interpolations/all.hpp
index 760e4f8fdf1..82c7d47acee 100644
--- a/ql/math/interpolations/all.hpp
+++ b/ql/math/interpolations/all.hpp
@@ -14,6 +14,7 @@
 #include <ql/math/interpolations/interpolation2d.hpp>
 #include <ql/math/interpolations/kernelinterpolation.hpp>
 #include <ql/math/interpolations/kernelinterpolation2d.hpp>
+#include <ql/math/interpolations/lagrangeinterpolation.hpp>
 #include <ql/math/interpolations/linearinterpolation.hpp>
 #include <ql/math/interpolations/loginterpolation.hpp>
 #include <ql/math/interpolations/mixedinterpolation.hpp>
diff --git a/ql/math/interpolations/lagrangeinterpolation.hpp b/ql/math/interpolations/lagrangeinterpolation.hpp
new file mode 100644
index 00000000000..913da668b73
--- /dev/null
+++ b/ql/math/interpolations/lagrangeinterpolation.hpp
@@ -0,0 +1,156 @@
+/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+
+/*
+ Copyright (C) 2016 Klaus Spanderen
+
+ This file is part of QuantLib, a free-software/open-source library
+ for financial quantitative analysts and developers - http://quantlib.org/
+
+ QuantLib is free software: you can redistribute it and/or modify it
+ under the terms of the QuantLib license.  You should have received a
+ copy of the license along with this program; if not, please email
+ <quantlib-dev@lists.sf.net>. The license is also available online at
+ <http://quantlib.org/license.shtml>.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE.  See the license for more details.
+*/
+
+#ifndef quantlib_lagrange_interpolation_hpp
+#define quantlib_lagrange_interpolation_hpp
+
+#include <ql/math/array.hpp>
+#include <ql/math/interpolation.hpp>
+
+#include <boost/make_shared.hpp>
+
+namespace QuantLib {
+    /*! References: J-P. Berrut and L.N. Trefethen,
+                    Barycentric Lagrange interpolation,
+                    SIAM Review, 46(3):501–517, 2004.
+        https://people.maths.ox.ac.uk/trefethen/barycentric.pdf
+    */
+
+    namespace detail {
+        class UpdatedYInterpolation {
+          public:
+            virtual ~UpdatedYInterpolation() {}
+            virtual Real value(const Array& yValues, Real x) const = 0;
+        };
+
+        template <class I1, class I2>
+        class LagrangeInterpolationImpl
+            : public Interpolation::templateImpl<I1,I2>,
+              public UpdatedYInterpolation {
+
+          public:
+            LagrangeInterpolationImpl(const I1& xBegin, const I1& xEnd,
+                                      const I2& yBegin)
+            : Interpolation::templateImpl<I1,I2>(xBegin, xEnd, yBegin),
+              n_(std::distance(xBegin, xEnd)),
+              lambda_(n_) {
+                #if defined(QL_EXTRA_SAFETY_CHECKS)
+                QL_REQUIRE(std::set<Real>(xBegin, xEnd).size() == n_,
+                        "x values must not contain duplicates");
+                #endif
+            }
+
+            void update() {
+                const Real cM1 = 4.0/(*(this->xEnd_-1) - *(this->xBegin_));
+
+                for (Size i=0; i < n_; ++i) {
+                    lambda_[i] = 1.0;
+
+                    const Real x_i = this->xBegin_[i];
+                    for (Size j=0; j < n_; ++j) {
+                        if (i != j)
+                            lambda_[i] *= cM1*(x_i-this->xBegin_[j]);
+                    }
+                    lambda_[i] = 1.0/lambda_[i];
+                }
+            }
+
+            Real value(Real x) const {
+                return _value(this->yBegin_, x);
+            }
+
+            Real derivative(Real x) const {
+                Real n=0.0, d=0.0, nd=0.0, dd=0.0;
+                for (Size i=0; i < n_; ++i) {
+                    const Real x_i = this->xBegin_[i];
+
+                    if (close_enough(x, x_i)) {
+                        Real p=0.0;
+                        for (Size j=0; j < n_; ++j)
+                            if (i != j) {
+                                p+=lambda_[j]/(x-this->xBegin_[j])
+                                    *(this->yBegin_[j] - this->yBegin_[i]);
+                            }
+                        return p/lambda_[i];
+                    }
+
+                    const Real alpha = lambda_[i]/(x-x_i);
+                    const Real alphad = -alpha/(x-x_i);
+                    n += alpha * this->yBegin_[i];
+                    d += alpha;
+                    nd += alphad * this->yBegin_[i];
+                    dd += alphad;
+                }
+                return (nd * d - n * dd)/(d*d);
+            }
+
+            Real primitive(Real) const {
+                QL_FAIL("LagrangeInterpolation primitive is not implemented");
+            }
+
+            Real secondDerivative(Real) const {
+                QL_FAIL("LagrangeInterpolation secondDerivative "
+                        "is not implemented");
+            }
+
+            Real value(const Array& y, Real x) const {
+                return _value(y.begin(), x);
+            }
+
+            private:
+              template <class Iterator>
+              Real _value(const Iterator& yBegin, Real x) const {
+                  Real n=0.0, d=0.0;
+                  for (Size i=0; i < n_; ++i) {
+                      const Real x_i = this->xBegin_[i];
+                      if (close_enough(x, x_i))
+                          return yBegin[i];
+
+                      const Real alpha = lambda_[i]/(x-x_i);
+                      n += alpha * yBegin[i];
+                      d += alpha;
+                  }
+                  return n/d;
+              }
+
+              const Size n_;
+              Array lambda_;
+        };
+    }
+
+    class LagrangeInterpolation : public Interpolation {
+      public:
+        template <class I1, class I2>
+        LagrangeInterpolation(const I1& xBegin, const I1& xEnd,
+                              const I2& yBegin) {
+            impl_ = boost::make_shared<detail::LagrangeInterpolationImpl<I1,I2> >(
+                xBegin, xEnd, yBegin);
+            impl_->update();
+        }
+
+        // interpolate with new set of y values for a new x value
+        Real value(const Array& y, Real x) const {
+            return boost::dynamic_pointer_cast<detail::UpdatedYInterpolation>
+                (impl_)->value(y, x);
+        }
+    };
+
+}
+
+#endif
diff --git a/ql/math/randomnumbers/Makefile.am b/ql/math/randomnumbers/Makefile.am
index de70b962aab..9c9a453334e 100644
--- a/ql/math/randomnumbers/Makefile.am
+++ b/ql/math/randomnumbers/Makefile.am
@@ -22,7 +22,8 @@ this_include_HEADERS = \
 	rngtraits.hpp \
 	seedgenerator.hpp \
 	sobolbrownianbridgersg.hpp \
-	sobolrsg.hpp
+	sobolrsg.hpp \
+	stochasticcollocationinvcdf.hpp
 
 libRandomNumbers_la_SOURCES = \
     faurersg.cpp \
@@ -35,7 +36,8 @@ libRandomNumbers_la_SOURCES = \
 	primitivepolynomials.cpp \
 	seedgenerator.cpp \
 	sobolbrownianbridgersg.cpp \
-	sobolrsg.cpp
+	sobolrsg.cpp \
+	stochasticcollocationinvcdf.cpp
 
 noinst_LTLIBRARIES = libRandomNumbers.la
 
diff --git a/ql/math/randomnumbers/all.hpp b/ql/math/randomnumbers/all.hpp
index e6371cff738..c830143a490 100644
--- a/ql/math/randomnumbers/all.hpp
+++ b/ql/math/randomnumbers/all.hpp
@@ -20,4 +20,5 @@
 #include <ql/math/randomnumbers/seedgenerator.hpp>
 #include <ql/math/randomnumbers/sobolbrownianbridgersg.hpp>
 #include <ql/math/randomnumbers/sobolrsg.hpp>
+#include <ql/math/randomnumbers/stochasticcollocationinvcdf.hpp>
 
diff --git a/ql/math/randomnumbers/stochasticcollocationinvcdf.cpp b/ql/math/randomnumbers/stochasticcollocationinvcdf.cpp
new file mode 100644
index 00000000000..0d1c8aae1cb
--- /dev/null
+++ b/ql/math/randomnumbers/stochasticcollocationinvcdf.cpp
@@ -0,0 +1,64 @@
+/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+
+/*
+
+ Copyright (C) 2016 Klaus Spanderen
+
+ This file is part of QuantLib, a free-software/open-source library
+ for financial quantitative analysts and developers - http://quantlib.org/
+
+ QuantLib is free software: you can redistribute it and/or modify it
+ under the terms of the QuantLib license.  You should have received a
+ copy of the license along with this program; if not, please email
+ <quantlib-dev@lists.sf.net>. The license is also available online at
+ <http://quantlib.org/license.shtml>.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE.  See the license for more details.
+*/
+
+/*! \file stochasticcollationcdf.cpp
+*/
+
+#include <ql/mathconstants.hpp>
+#include <ql/math/integrals/gaussianquadratures.hpp>
+#include <ql/math/randomnumbers/stochasticcollocationinvcdf.hpp>
+
+namespace QuantLib {
+
+    namespace {
+        Disposable<Array> g(Real sigma, const Array& x,
+                            const boost::function<Real(Real)>& invCDF) {
+
+            Array y(x.size());
+            const CumulativeNormalDistribution normalCDF;
+
+            for (Size i=0, n=x.size(); i < n; ++i) {
+                y[i] = invCDF(normalCDF(x[i]/sigma));
+            }
+
+            return y;
+        }
+    }
+
+    StochasticCollocationInvCDF::StochasticCollocationInvCDF(
+        const boost::function<Real(Real)>& invCDF,
+        Size lagrangeOrder, Real pMax, Real pMin)
+    : x_(M_SQRT2*GaussHermiteIntegration(lagrangeOrder).x()),
+      sigma_( (pMax != Null<Real>())
+              ? x_.back() / InverseCumulativeNormal()(pMax)
+              : (pMin != Null<Real>())
+                  ? x_.front() / InverseCumulativeNormal()(pMin)
+                  : 1.0),
+      y_(g(sigma_, x_, invCDF)),
+      interpl_(x_.begin(), x_.end(), y_.begin()) {
+    }
+
+    Real StochasticCollocationInvCDF::value(Real x) const {
+        return interpl_(x*sigma_, true);
+    }
+    Real StochasticCollocationInvCDF::operator()(Real u) const {
+        return value(InverseCumulativeNormal()(u));
+    }
+}
diff --git a/ql/math/randomnumbers/stochasticcollocationinvcdf.hpp b/ql/math/randomnumbers/stochasticcollocationinvcdf.hpp
new file mode 100755
index 00000000000..b7aec342428
--- /dev/null
+++ b/ql/math/randomnumbers/stochasticcollocationinvcdf.hpp
@@ -0,0 +1,63 @@
+/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+
+/*
+
+ Copyright (C) 2016 Klaus Spanderen
+
+ This file is part of QuantLib, a free-software/open-source library
+ for financial quantitative analysts and developers - http://quantlib.org/
+
+ QuantLib is free software: you can redistribute it and/or modify it
+ under the terms of the QuantLib license.  You should have received a
+ copy of the license along with this program; if not, please email
+ <quantlib-dev@lists.sf.net>. The license is also available online at
+ <http://quantlib.org/license.shtml>.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE.  See the license for more details.
+*/
+
+/*! \file stochasticcollationcdf.hpp
+    Stochastic collocation inverse cumulative distribution function
+*/
+
+#ifndef quantlib_stochastic_collation_inv_cdf_hpp
+#define quantlib_stochastic_collation_inv_cdf_hpp
+
+#include <ql/math/distributions/normaldistribution.hpp>
+#include <ql/math/interpolations/lagrangeinterpolation.hpp>
+
+#include <boost/function.hpp>
+#include <functional>
+
+namespace QuantLib {
+    //! Stochastic collocation inverse cumulative distribution function
+
+    /*! References:
+        L.A. Grzelak, J.A.S. Witteveen, M.Suárez-Taboada, C.W. Oosterlee,
+        The Stochastic Collocation Monte Carlo Sampler: Highly efficient
+        sampling from “expensive” distributions
+        http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2529691
+     */
+
+    class StochasticCollocationInvCDF : public std::unary_function<Real,Real> {
+      public:
+        StochasticCollocationInvCDF(
+            const boost::function<Real(Real)>& invCDF,
+            Size lagrangeOrder,
+            Real pMax = Null<Real>(),
+            Real pMin = Null<Real>());
+
+        Real value(Real x) const;
+        Real operator()(Real u) const;
+
+      private:
+        const Array x_;
+        const Volatility sigma_;
+        const Array y_;
+        const LagrangeInterpolation interpl_;
+    };
+}
+
+#endif
diff --git a/ql/methods/finitedifferences/operators/Makefile.am b/ql/methods/finitedifferences/operators/Makefile.am
index e134893e8cb..7782c43de5a 100644
--- a/ql/methods/finitedifferences/operators/Makefile.am
+++ b/ql/methods/finitedifferences/operators/Makefile.am
@@ -12,6 +12,7 @@ this_include_HEADERS = \
 	fdmhestonop.hpp \
 	fdmhullwhiteop.hpp \
 	fdmlinearopcomposite.hpp \
+	fdmornsteinuhlenbeckop.hpp \
 	fdmlinearop.hpp \
 	fdmlinearopiterator.hpp \
 	fdmlinearoplayout.hpp \
@@ -30,6 +31,7 @@ libFdmOperators_la_SOURCES = \
 	fdmhestonop.cpp \
 	fdmhullwhiteop.cpp \
 	fdmlinearoplayout.cpp \
+	fdmornsteinuhlenbeckop.cpp \
 	firstderivativeop.cpp \
 	ninepointlinearop.cpp \
 	secondderivativeop.cpp \
diff --git a/ql/methods/finitedifferences/operators/all.hpp b/ql/methods/finitedifferences/operators/all.hpp
index 918c467285d..6177b9c31ea 100644
--- a/ql/methods/finitedifferences/operators/all.hpp
+++ b/ql/methods/finitedifferences/operators/all.hpp
@@ -9,6 +9,7 @@
 #include <ql/methods/finitedifferences/operators/fdmhestonop.hpp>
 #include <ql/methods/finitedifferences/operators/fdmhullwhiteop.hpp>
 #include <ql/methods/finitedifferences/operators/fdmlinearopcomposite.hpp>
+#include <ql/methods/finitedifferences/operators/fdmornsteinuhlenbeckop.hpp>
 #include <ql/methods/finitedifferences/operators/fdmlinearop.hpp>
 #include <ql/methods/finitedifferences/operators/fdmlinearopiterator.hpp>
 #include <ql/methods/finitedifferences/operators/fdmlinearoplayout.hpp>
diff --git a/ql/methods/finitedifferences/operators/fdmornsteinuhlenbeckop.cpp b/ql/methods/finitedifferences/operators/fdmornsteinuhlenbeckop.cpp
new file mode 100644
index 00000000000..03d3298698b
--- /dev/null
+++ b/ql/methods/finitedifferences/operators/fdmornsteinuhlenbeckop.cpp
@@ -0,0 +1,122 @@
+/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+
+/*
+ Copyright (C) 2016 Klaus Spanderen
+
+ This file is part of QuantLib, a free-software/open-source library
+ for financial quantitative analysts and developers - http://quantlib.org/
+
+ QuantLib is free software: you can redistribute it and/or modify it
+ under the terms of the QuantLib license.  You should have received a
+ copy of the license along with this program; if not, please email
+ <quantlib-dev@lists.sf.net>. The license is also available online at
+ <http://quantlib.org/license.shtml>.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE.  See the license for more details.
+*/
+
+/*! \file fdmornsteinuhlenbeckop.cpp
+*/
+
+#include <ql/math/functional.hpp>
+#include <ql/processes/ornsteinuhlenbeckprocess.hpp>
+#include <ql/termstructures/yieldtermstructure.hpp>
+#include <ql/methods/finitedifferences/meshers/fdmmesher.hpp>
+#include <ql/methods/finitedifferences/operators/fdmlinearoplayout.hpp>
+#include <ql/methods/finitedifferences/operators/fdmornsteinuhlenbeckop.hpp>
+#include <ql/methods/finitedifferences/operators/secondderivativeop.hpp>
+
+namespace QuantLib {
+
+    FdmOrnsteinUhlenbackOp::FdmOrnsteinUhlenbackOp(
+            const boost::shared_ptr<FdmMesher>& mesher,
+            const boost::shared_ptr<OrnsteinUhlenbeckProcess>& process,
+            const boost::shared_ptr<YieldTermStructure>& rTS,
+            const FdmBoundaryConditionSet& bcSet,
+            Size direction)
+    : mesher_   (mesher),
+      process_  (process),
+      rTS_      (rTS),
+      bcSet_    (bcSet),
+      direction_(direction),
+      m_        (direction, mesher),
+      mapX_     (direction, mesher)  {
+
+        const boost::shared_ptr<FdmLinearOpLayout> layout=mesher_->layout();
+
+        Array drift(layout->size());
+        const Array x(mesher_->locations(direction));
+
+        for (FdmLinearOpIterator iter=layout->begin(), endIter=layout->end();
+             iter!=endIter; ++iter) {
+            const Size i = iter.index();
+            drift[i] = process_->drift(0.0, x[i]);
+        }
+
+        m_.axpyb(drift, FirstDerivativeOp(direction, mesher),
+            SecondDerivativeOp(direction, mesher)
+                .mult(0.5*square<Real>()(process_->volatility())
+                      *Array(mesher->layout()->size(), 1.0)), Array());
+    }
+
+    Size FdmOrnsteinUhlenbackOp::size() const {
+        return mesher_->layout()->dim().size();;
+    }
+
+    void FdmOrnsteinUhlenbackOp::setTime(Time t1, Time t2) {
+        const Rate r = rTS_->forwardRate(t1, t2, Continuous).rate();
+
+        mapX_.axpyb(Array(), m_, m_, Array(1, -r));
+    }
+
+    Disposable<Array> FdmOrnsteinUhlenbackOp::apply(const Array& r) const {
+        return mapX_.apply(r);
+    }
+
+    Disposable<Array> FdmOrnsteinUhlenbackOp::apply_mixed(
+        const Array& r) const {
+
+        Array retVal(r.size(), 0.0);
+        return retVal;
+    }
+
+    Disposable<Array> FdmOrnsteinUhlenbackOp::apply_direction(
+        Size direction, const Array& r) const {
+
+        if (direction == direction_) {
+            return mapX_.apply(r);
+        }
+        else {
+            Array retVal(r.size(), 0.0);
+            return retVal;
+        }
+    }
+
+    Disposable<Array> FdmOrnsteinUhlenbackOp::solve_splitting(
+        Size direction, const Array& r, Real a) const {
+
+        if (direction == direction_) {
+            return mapX_.solve_splitting(r, a, 1.0);
+        }
+        else {
+            Array retVal(r);
+            return retVal;
+        }
+    }
+
+    Disposable<Array> FdmOrnsteinUhlenbackOp::preconditioner(
+        const Array& r, Real dt) const {
+        return solve_splitting(direction_, r, dt);
+    }
+
+#if !defined(QL_NO_UBLAS_SUPPORT)
+    Disposable<std::vector<SparseMatrix> >
+    FdmOrnsteinUhlenbackOp::toMatrixDecomp() const {
+        std::vector<SparseMatrix> retVal(1, mapX_.toMatrix());
+        return retVal;
+    }
+#endif
+
+}
diff --git a/ql/methods/finitedifferences/operators/fdmornsteinuhlenbeckop.hpp b/ql/methods/finitedifferences/operators/fdmornsteinuhlenbeckop.hpp
new file mode 100644
index 00000000000..17aa7ff859e
--- /dev/null
+++ b/ql/methods/finitedifferences/operators/fdmornsteinuhlenbeckop.hpp
@@ -0,0 +1,73 @@
+/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+
+/*
+ Copyright (C) 2016 Klaus Spanderen
+
+ This file is part of QuantLib, a free-software/open-source library
+ for financial quantitative analysts and developers - http://quantlib.org/
+
+ QuantLib is free software: you can redistribute it and/or modify it
+ under the terms of the QuantLib license.  You should have received a
+ copy of the license along with this program; if not, please email
+ <quantlib-dev@lists.sf.net>. The license is also available online at
+ <http://quantlib.org/license.shtml>.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE.  See the license for more details.
+*/
+
+/*! \file fdmornsteinuhlenbeckop.hpp
+    \brief Ornstein Uhlenbeck process
+*/
+
+#ifndef quantlib_fdm_ornstein_uhlenback_op_hpp
+#define quantlib_fdm_ornstein_uhlenback_op_hpp
+
+#include <ql/methods/finitedifferences/operators/firstderivativeop.hpp>
+#include <ql/methods/finitedifferences/operators/triplebandlinearop.hpp>
+#include <ql/methods/finitedifferences/operators/fdmlinearopcomposite.hpp>
+#include <ql/methods/finitedifferences/utilities/fdmboundaryconditionset.hpp>
+
+
+namespace QuantLib {
+
+    class FdmMesher;
+    class YieldTermStructure;
+    class OrnsteinUhlenbeckProcess;
+
+    class FdmOrnsteinUhlenbackOp : public FdmLinearOpComposite {
+      public:
+        FdmOrnsteinUhlenbackOp(
+            const boost::shared_ptr<FdmMesher>& mesher,
+            const boost::shared_ptr<OrnsteinUhlenbeckProcess>& p,
+            const boost::shared_ptr<YieldTermStructure>& rTS,
+            const FdmBoundaryConditionSet& bcSet,
+            Size direction = 0);
+
+        Size size() const;
+        void setTime(Time t1, Time t2);
+
+        Disposable<Array> apply(const Array& r) const;
+        Disposable<Array> apply_mixed(const Array& r) const;
+
+        Disposable<Array> apply_direction(Size direction,
+                                          const Array& r) const;
+        Disposable<Array> solve_splitting(Size direction,
+                                          const Array& r, Real s) const;
+        Disposable<Array> preconditioner(const Array& r, Real s) const;
+
+#if !defined(QL_NO_UBLAS_SUPPORT)
+        Disposable<std::vector<SparseMatrix> > toMatrixDecomp() const;
+#endif
+      private:
+        const boost::shared_ptr<FdmMesher> mesher_;
+        const boost::shared_ptr<OrnsteinUhlenbeckProcess> process_;
+        const boost::shared_ptr<YieldTermStructure> rTS_;
+        const FdmBoundaryConditionSet bcSet_;
+        const Size direction_;
+
+        TripleBandLinearOp m_, mapX_;
+    };
+}
+#endif
diff --git a/ql/methods/finitedifferences/solvers/fdm2dblackscholessolver.cpp b/ql/methods/finitedifferences/solvers/fdm2dblackscholessolver.cpp
index 78aaca5f213..a12f5fd377a 100644
--- a/ql/methods/finitedifferences/solvers/fdm2dblackscholessolver.cpp
+++ b/ql/methods/finitedifferences/solvers/fdm2dblackscholessolver.cpp
@@ -107,4 +107,13 @@ namespace QuantLib {
         return (solver_->derivativeYY(x, y)
                 -solver_->derivativeY(x, y))/(v*v);
     }
+
+    Real Fdm2dBlackScholesSolver::gammaXYat(Real u, Real v) const {
+        calculate();
+
+        const Real x = std::log(u);
+        const Real y = std::log(v);
+
+        return solver_->derivativeXY(x, y)/(u*v);
+    }
 }
diff --git a/ql/methods/finitedifferences/solvers/fdm2dblackscholessolver.hpp b/ql/methods/finitedifferences/solvers/fdm2dblackscholessolver.hpp
index 41d3afd3fa8..3c23203f0aa 100644
--- a/ql/methods/finitedifferences/solvers/fdm2dblackscholessolver.hpp
+++ b/ql/methods/finitedifferences/solvers/fdm2dblackscholessolver.hpp
@@ -49,6 +49,7 @@ namespace QuantLib {
         Real deltaYat(Real x, Real y) const;
         Real gammaXat(Real x, Real y) const;
         Real gammaYat(Real x, Real y) const;
+        Real gammaXYat(Real x, Real y) const;
 
       protected:
         void performCalculations() const;
diff --git a/ql/pricingengines/basket/fd2dblackscholesvanillaengine.cpp b/ql/pricingengines/basket/fd2dblackscholesvanillaengine.cpp
index a2110afe6f4..33bf64dcc02 100644
--- a/ql/pricingengines/basket/fd2dblackscholesvanillaengine.cpp
+++ b/ql/pricingengines/basket/fd2dblackscholesvanillaengine.cpp
@@ -41,6 +41,8 @@ namespace QuantLib {
       xGrid_(xGrid), yGrid_(yGrid), tGrid_(tGrid),
       dampingSteps_(dampingSteps),
       schemeDesc_(schemeDesc) {
+        registerWith(p1);
+        registerWith(p2);
     }
 
     void Fd2dBlackScholesVanillaEngine::calculate() const {
@@ -95,6 +97,9 @@ namespace QuantLib {
         const Real y = p2_->x0();
 
         results_.value = solver->valueAt(x, y);
+        results_.delta = solver->deltaXat(x, y) + solver->deltaYat(x, y);
+        results_.gamma = solver->gammaXat(x, y) + solver->gammaYat(x, y)
+             + 2*solver->gammaXYat(x, y);
         results_.theta = solver->thetaAt(x, y);
     }
 }
diff --git a/ql/processes/ornsteinuhlenbeckprocess.cpp b/ql/processes/ornsteinuhlenbeckprocess.cpp
index e7d4a189359..eddb653054b 100644
--- a/ql/processes/ornsteinuhlenbeckprocess.cpp
+++ b/ql/processes/ornsteinuhlenbeckprocess.cpp
@@ -28,12 +28,11 @@ namespace QuantLib {
                                                        Real x0,
                                                        Real level)
     : x0_(x0), speed_(speed), level_(level), volatility_(vol) {
-        QL_REQUIRE(speed_ >= 0.0, "negative speed given");
         QL_REQUIRE(volatility_ >= 0.0, "negative volatility given");
     }
 
     Real OrnsteinUhlenbeckProcess::variance(Time, Real, Time dt) const {
-        if (speed_ < std::sqrt(QL_EPSILON)) {
+        if (std::fabs(speed_) < std::sqrt(QL_EPSILON)) {
              // algebraic limit for small speed
             return volatility_*volatility_*dt;
         } else {
diff --git a/test-suite/Makefile.am b/test-suite/Makefile.am
index 8339628262d..e5863661b3c 100644
--- a/test-suite/Makefile.am
+++ b/test-suite/Makefile.am
@@ -92,6 +92,7 @@ QL_TESTS = \
 	mersennetwister.hpp mersennetwister.cpp \
 	money.hpp money.cpp \
 	noarbsabr.hpp noarbsabr.cpp \
+	normalclvmodel.hpp normalclvmodel.cpp \
 	nthtodefault.hpp nthtodefault.cpp \
 	numericaldifferentiation.hpp numericaldifferentiation.cpp \
 	observable.hpp observable.cpp \
diff --git a/test-suite/basketoption.cpp b/test-suite/basketoption.cpp
index 59cceefe1f4..821846b2792 100644
--- a/test-suite/basketoption.cpp
+++ b/test-suite/basketoption.cpp
@@ -33,7 +33,9 @@
 #include <ql/termstructures/volatility/equityfx/blackconstantvol.hpp>
 #include <ql/pricingengines/basket/fd2dblackscholesvanillaengine.hpp>
 #include <ql/utilities/dataformatters.hpp>
+
 #include <boost/progress.hpp>
+#include <boost/make_shared.hpp>
 #include <boost/bind.hpp>
 
 using namespace QuantLib;
@@ -946,8 +948,84 @@ void BasketOptionTest::testOddSamples() {
     }
 }
 
+void BasketOptionTest::test2DPDEGreeks() {
+
+    BOOST_TEST_MESSAGE("Testing Greeks of 2dim PDE engine...");
+
+    const Real s1 = 100;
+    const Real s2 = 100;
+    const Real r = 0.013;
+    const Volatility v = 0.2;
+    const Real rho = 0.5;
+    const Real strike = s1-s2;
+    const Size maturityInDays = 1095;
+
+    const DayCounter dc = Actual365Fixed();
+    const Date today = Date::todaysDate();
+    const Date maturity = today + maturityInDays;
+
+    const boost::shared_ptr<SimpleQuote> spot1(
+        boost::make_shared<SimpleQuote>(s1));
+    const boost::shared_ptr<SimpleQuote> spot2(
+        boost::make_shared<SimpleQuote>(s2));
+
+    const Handle<YieldTermStructure> rTS(flatRate(today, r, dc));
+    const Handle<BlackVolTermStructure> vTS(flatVol(today, v, dc));
+
+    const boost::shared_ptr<BlackProcess> p1(
+        boost::make_shared<BlackProcess>(Handle<Quote>(spot1), rTS, vTS));
+
+    const boost::shared_ptr<BlackProcess> p2(
+        boost::make_shared<BlackProcess>(Handle<Quote>(spot2), rTS, vTS));
+
+    BasketOption option(
+        boost::make_shared<SpreadBasketPayoff>(
+            boost::make_shared<PlainVanillaPayoff>(Option::Call, strike)),
+        boost::make_shared<EuropeanExercise>(maturity));
+
+    option.setPricingEngine(
+        boost::make_shared<Fd2dBlackScholesVanillaEngine>(p1, p2, rho));
+
+    const Real calculatedDelta = option.delta();
+    const Real calculatedGamma = option.gamma();
+
+    option.setPricingEngine(boost::make_shared<KirkEngine>(p1, p2, rho));
+
+    const Real eps = 1.0;
+    const Real npv = option.NPV();
+
+    spot1->setValue(s1 + eps);
+    spot2->setValue(s2 + eps);
+    const Real npvUp = option.NPV();
+
+    spot1->setValue(s1 - eps);
+    spot2->setValue(s2 - eps);
+    const Real npvDown = option.NPV();
+
+    const Real expectedDelta = (npvUp - npvDown)/(2*eps);
+    const Real expectedGamma = (npvUp + npvDown - 2*npv)/(eps*eps);
+
+    const Real tol = 0.0005;
+    if (std::fabs(expectedDelta - calculatedDelta) > tol) {
+        BOOST_FAIL("failed to reproduce delta with 2dim PDE"
+                   << QL_FIXED << std::setprecision(8)
+                   << "\n    calculated: " << calculatedDelta
+                   << "\n    expected:   " << expectedDelta
+                   << "\n    tolerance:  " << tol);
+    }
+
+    if (std::fabs(expectedGamma - calculatedGamma) > tol) {
+        BOOST_FAIL("failed to reproduce delta with 2dim PDE"
+                   << QL_FIXED << std::setprecision(8)
+                   << "\n    calculated: " << calculatedGamma
+                   << "\n    expected:   " << expectedGamma
+                   << "\n    tolerance:  " << tol);
+    }
+}
+
 test_suite* BasketOptionTest::suite() {
     test_suite* suite = BOOST_TEST_SUITE("Basket option tests");
+
     suite->add(QUANTLIB_TEST_CASE(&BasketOptionTest::testEuroTwoValues));
     // FLOATING_POINT_EXCEPTION
     suite->add(QUANTLIB_TEST_CASE(
@@ -962,6 +1040,7 @@ test_suite* BasketOptionTest::suite() {
                 ((i+1)*LENGTH(oneDataValues))/nTestCases)));
     }
     suite->add(QUANTLIB_TEST_CASE(&BasketOptionTest::testOddSamples));
+    suite->add(QUANTLIB_TEST_CASE(&BasketOptionTest::test2DPDEGreeks));
 
     return suite;
 }
diff --git a/test-suite/basketoption.hpp b/test-suite/basketoption.hpp
index ad486f83f39..8d9916eaf10 100644
--- a/test-suite/basketoption.hpp
+++ b/test-suite/basketoption.hpp
@@ -32,6 +32,7 @@ class BasketOptionTest {
     static void testTavellaValues();
     static void testOneDAmericanValues(unsigned from, unsigned to);
     static void testOddSamples();
+    static void test2DPDEGreeks();
     static boost::unit_test_framework::test_suite* suite();
 };
 
diff --git a/test-suite/distributions.cpp b/test-suite/distributions.cpp
index a65d6569551..771573e81ad 100644
--- a/test-suite/distributions.cpp
+++ b/test-suite/distributions.cpp
@@ -5,6 +5,8 @@
  Copyright (C) 2003 StatPro Italia srl
  Copyright (C) 2005 Gary Kennedy
  Copyright (C) 2013 Fabien Le Floc'h
+ Copyright (C) 2016 Klaus Spanderen
+
 
  This file is part of QuantLib, a free-software/open-source library
  for financial quantitative analysts and developers - http://quantlib.org/
@@ -25,10 +27,14 @@
 #include <ql/math/distributions/normaldistribution.hpp>
 #include <ql/math/distributions/bivariatenormaldistribution.hpp>
 #include <ql/math/distributions/bivariatestudenttdistribution.hpp>
+#include <ql/math/distributions/chisquaredistribution.hpp>
 #include <ql/math/distributions/poissondistribution.hpp>
+#include <ql/math/randomnumbers/stochasticcollocationinvcdf.hpp>
 #include <ql/math/comparison.hpp>
 #include <ql/math/functional.hpp>
 
+#include <boost/math/distributions/non_central_chi_squared.hpp>
+
 using namespace QuantLib;
 using namespace boost::unit_test_framework;
 
@@ -624,6 +630,86 @@ void DistributionTest::testBivariateCumulativeStudentVsBivariate() {
     }
 }
     
+
+namespace {
+    class InverseNonCentralChiSquared : public std::unary_function<Real,Real> {
+      public:
+        InverseNonCentralChiSquared(Real df, Real ncp)
+        : dist_(df, ncp) {}
+
+        Real operator()(Real x) const {
+            return boost::math::quantile(dist_, x);
+        }
+      private:
+        const boost::math::non_central_chi_squared_distribution<Real> dist_;
+    };
+}
+
+void DistributionTest::testInvCDFviaStochasticCollocation() {
+    BOOST_TEST_MESSAGE(
+        "Testing inverse CDF based on testInvCDFviaStochasticCollocation...");
+
+    const Real k = 3.0;
+    const Real lambda = 1.0;
+
+    const InverseCumulativeNormal invNormalCDF;
+    const CumulativeNormalDistribution normalCDF;
+    const InverseNonCentralChiSquared invCDF(k, lambda);
+
+    const StochasticCollocationInvCDF scInvCDF10(invCDF, 10);
+
+    // low precision
+    for (Real x=-3.0; x < 3.0; x+=0.1) {
+        const Real u = normalCDF(x);
+
+        const Real calculated1 = scInvCDF10(u);
+        const Real calculated2 = scInvCDF10.value(x);
+        const Real expected = invCDF(u);
+
+        if (std::fabs(calculated1 - calculated2) > 1e-6) {
+            BOOST_FAIL("Failed to reproduce equal stochastic collocation "
+                       "inverse CDF" <<
+                       "\n    x: " << x <<
+                       "\n    calculated via normal distribution : "
+                           << calculated2 <<
+                       "\n    calculated via uniform distribution: "
+                           << calculated1 <<
+                       "\n    diff: " << calculated1 - calculated2);
+        }
+
+        const Real tol = 1e-2;
+        if (std::fabs(calculated2 - expected) > tol) {
+            BOOST_FAIL("Failed to reproduce invCDF with "
+                       "stochastic collocation method" <<
+                       "\n    x: " << x <<
+                       "\n    invCDF  :" << expected <<
+                       "\n    scInvCDF: " << calculated2 <<
+                       "\n    diff    : " << std::fabs(expected-calculated2) <<
+                       "\n    tol     : " << tol);
+        }
+    }
+
+    // high precision
+    const StochasticCollocationInvCDF scInvCDF30(invCDF, 30, 0.9999999);
+    for (Real x=-4.0; x < 4.0; x+=0.1) {
+        const Real u = normalCDF(x);
+
+        const Real expected = invCDF(u);
+        const Real calculated = scInvCDF30(u);
+
+        const Real tol = 1e-6;
+        if (std::fabs(calculated - expected) > tol) {
+            BOOST_FAIL("Failed to reproduce invCDF with "
+                       "stochastic collocation method" <<
+                       "\n    x: " << x <<
+                       "\n    invCDF  :" << expected <<
+                       "\n    scInvCDF: " << calculated <<
+                       "\n    diff    : " << std::fabs(expected-calculated) <<
+                       "\n    tol     : " << tol);
+        }
+    }
+}
+
 test_suite* DistributionTest::suite() {
     test_suite* suite = BOOST_TEST_SUITE("Distribution tests");
     suite->add(QUANTLIB_TEST_CASE(&DistributionTest::testNormal));
@@ -636,6 +722,9 @@ test_suite* DistributionTest::suite() {
                           &DistributionTest::testBivariateCumulativeStudent));
     suite->add(QUANTLIB_TEST_CASE(
                &DistributionTest::testBivariateCumulativeStudentVsBivariate));
+    suite->add(QUANTLIB_TEST_CASE(
+                   &DistributionTest::testInvCDFviaStochasticCollocation));
+
     return suite;
 }
 
diff --git a/test-suite/distributions.hpp b/test-suite/distributions.hpp
index d571a56bcd0..6b100492df7 100644
--- a/test-suite/distributions.hpp
+++ b/test-suite/distributions.hpp
@@ -3,6 +3,7 @@
 /*
  Copyright (C) 2003 Ferdinando Ametrano
  Copyright (C) 2003 RiskMap srl
+ Copyright (C) 2016 Klaus Spanderen
 
  This file is part of QuantLib, a free-software/open-source library
  for financial quantitative analysts and developers - http://quantlib.org/
@@ -35,6 +36,7 @@ class DistributionTest {
     static void testInverseCumulativePoisson();
     static void testBivariateCumulativeStudent();
     static void testBivariateCumulativeStudentVsBivariate();
+    static void testInvCDFviaStochasticCollocation();
     static boost::unit_test_framework::test_suite* suite();
 };
 
diff --git a/test-suite/interpolations.cpp b/test-suite/interpolations.cpp
index 7094b9d2d51..6e62df7ba5e 100644
--- a/test-suite/interpolations.cpp
+++ b/test-suite/interpolations.cpp
@@ -26,6 +26,7 @@
 #include <ql/utilities/dataformatters.hpp>
 #include <ql/utilities/null.hpp>
 #include <ql/math/interpolations/linearinterpolation.hpp>
+#include <ql/math/interpolations/bicubicsplineinterpolation.hpp>
 #include <ql/math/interpolations/backwardflatinterpolation.hpp>
 #include <ql/math/interpolations/forwardflatinterpolation.hpp>
 #include <ql/math/interpolations/cubicinterpolation.hpp>
@@ -33,7 +34,7 @@
 #include <ql/math/interpolations/sabrinterpolation.hpp>
 #include <ql/math/interpolations/kernelinterpolation.hpp>
 #include <ql/math/interpolations/kernelinterpolation2d.hpp>
-#include <ql/math/interpolations/bicubicsplineinterpolation.hpp>
+#include <ql/math/interpolations/lagrangeinterpolation.hpp>
 #include <ql/math/integrals/simpsonintegral.hpp>
 #include <ql/math/kernelfunctions.hpp>
 #include <ql/math/functional.hpp>
@@ -43,7 +44,7 @@
 #include <ql/experimental/volatility/noarbsabrinterpolation.hpp>
 #include <boost/foreach.hpp>
 #include <boost/assign/std/vector.hpp>
-
+#include <boost/math/special_functions/fpclassify.hpp>
 
 using namespace QuantLib;
 using namespace boost::unit_test_framework;
@@ -1981,6 +1982,177 @@ void InterpolationTest::testTransformations() {
 
 }
 
+namespace {
+    Real lagrangeTestFct(Real x) {
+        return std::fabs(x) + 0.5*x - x*x;
+    }
+}
+
+void InterpolationTest::testLagrangeInterpolation() {
+
+    BOOST_TEST_MESSAGE("Testing Lagrange Interpolation...");
+
+    const Real x[] = {-1.0 , -0.5, -0.25, 0.1, 0.4, 0.75, 0.96};
+    Array y(LENGTH(x));
+    std::transform(x, x+LENGTH(x), y.begin(), &lagrangeTestFct);
+
+    LagrangeInterpolation interpl(&x[0], x+LENGTH(x), y.begin());
+
+    // reference results are taken from R package pracma
+    const Real references[] = {
+        -0.5000000000000000,-0.5392414024347419,-0.5591485962711904,
+        -0.5629199661387594,-0.5534414777017116,-0.5333043347921566,
+        -0.5048221831582063,-0.4700478608272949,-0.4307896950846587,
+        -0.3886273460669714,-0.3449271969711449,-0.3008572908782903,
+        -0.2574018141928359,-0.2153751266968088,-0.1754353382192734,
+        -0.1380974319209344,-0.1037459341938971,-0.0726471311765894,
+        -0.0449608318838433,-0.0207516779521373,0.0000000000000000,
+        0.0173877793964286,0.0315691961126723,0.0427562482700356,
+        0.0512063534145595,0.0572137590808174,0.0611014067405497,
+        0.0632132491361394,0.0639070209989264,0.0635474631523613,
+        0.0625000000000000,0.0611248703983366,0.0597717119144768,
+        0.0587745984686508,0.0584475313615655,0.0590803836865967,
+        0.0609352981268212,0.0642435381368876,0.0692027925097279,
+        0.0759749333281079,0.0846842273010179,0.0954160004849021,
+        0.1082157563897290,0.1230887474699003,0.1400000000000001,
+        0.1588747923353829,0.1795995865576031,0.2020234135046815,
+        0.2259597111862140,0.2511886165833182,0.2774597108334206,
+        0.3044952177998833,0.3319936560264689,0.3596339440766487,
+        0.3870799592577457,0.4139855497299214,0.4400000000000001,
+        0.4647739498001331,0.4879657663513030,0.5092483700116673,
+        0.5283165133097421,0.5448945133624253,0.5587444376778583,
+        0.5696747433431296,0.5775493695968156,0.5822972837863635,
+        0.5839224807103117,0.5825144353453510,0.5782590089582251,
+        0.5714498086024714,0.5625000000000000,0.5519545738075141,
+        0.5405030652677689,0.5289927272456703,0.5184421566492137,
+        0.5100553742352614,0.5052363578001620,0.5056040287552059,
+        0.5130076920869246
+    };
+
+    const Real tol = 50*QL_EPSILON;
+    for (Size i=0; i < 79; ++i) {
+        const Real xx = -1.0 + i*0.025;
+        const Real calculated = interpl(xx);
+        if (   boost::math::isnan(calculated)
+            || std::fabs(references[i] - calculated) > tol) {
+            BOOST_FAIL("failed to reproduce the Lagrange interplation"
+                    << "\n    x         : " << xx
+                    << "\n    calculated: " << calculated
+                    << "\n    expected  : " << references[i]);
+        }
+    }
+}
+
+void InterpolationTest::testLagrangeInterpolationAtSupportPoint() {
+    BOOST_TEST_MESSAGE(
+        "Testing Lagrange Interpolation at supporting points...");
+
+    const Size n=5;
+    Array x(n), y(n);
+    for (Size i=0; i < n; ++i) {
+        x[i] = i/Real(n);
+        y[i] = 1.0/(1.0 - x[i]);
+    }
+    LagrangeInterpolation interpl(x.begin(), x.end(), y.begin());
+
+    const Real relTol = 5e-12;
+
+    for (Size i=1; i < n-1; ++i) {
+        for (Real z = x[i] - 100*QL_EPSILON;
+            z < x[i] + 100*QL_EPSILON; z+=2*QL_EPSILON) {
+            const Real expected = 1.0/(1.0 - x[i]);
+            const Real calculated = interpl(z);
+
+            if (   boost::math::isnan(calculated)
+                || std::fabs(expected - calculated) > relTol) {
+                BOOST_FAIL("failed to reproduce the Lagrange interplation"
+                        << "\n    x         : " << z
+                        << "\n    calculated: " << calculated
+                        << "\n    expected  : " << expected);
+            }
+        }
+    }
+}
+
+void InterpolationTest::testLagrangeInterpolationDerivative() {
+    BOOST_TEST_MESSAGE(
+        "Testing Lagrange Interpolation derivatives...");
+
+    Array x(5), y(5);
+    x[0] = -1.0; y[0] = 2.0;
+    x[1] = -0.3; y[1] = 3.0;
+    x[2] =  0.1; y[2] = 6.0;
+    x[3] =  0.3; y[3] = 3.0;
+    x[4] =  0.9; y[4] =-1.0;
+
+    LagrangeInterpolation interpl(x.begin(), x.end(), y.begin());
+
+    const Real eps = std::sqrt(QL_EPSILON);
+    for (Real x=-1.0; x <= 0.9; x+=0.01) {
+        const Real calculated = interpl.derivative(x, true);
+        const Real expected = (interpl(x+eps, true)
+            - interpl(x-eps, true))/(2*eps);
+
+        if (   boost::math::isnan(calculated)
+            || std::fabs(expected - calculated) > 25*eps) {
+            BOOST_FAIL("failed to reproduce the Lagrange"
+                    " interplation derivative"
+                    << "\n    x         : " << x
+                    << "\n    calculated: " << calculated
+                    << "\n    expected  : " << expected);
+        }
+    }
+}
+
+void InterpolationTest::testLagrangeInterpolationOnChebyshevPoints() {
+    BOOST_TEST_MESSAGE(
+        "Testing Lagrange Interpolation on Chebyshev points...");
+
+    // Test example taken from
+    // J.P. Berrut, L.N. Trefethen, Barycentric Lagrange Interpolation
+    // https://people.maths.ox.ac.uk/trefethen/barycentric.pdf
+
+    const Size n=50;
+    Array x(n+1), y(n+1);
+    for (Size i=0; i <= n; ++i) {
+        // Chebyshev points
+        x[i] = std::cos( (2*i+1)*M_PI/(2*n+2) );
+        y[i] = std::exp(x[i])/std::cos(x[i]);
+    }
+
+    LagrangeInterpolation interpl(x.begin(), x.end(), y.begin());
+
+    const Real tol = 1e-13;
+
+    for (Real x=-1.0; x <= 1.0; x+=0.01) {
+        const Real calculated = interpl(x, true);
+        const Real expected = std::exp(x)/std::cos(x);
+
+        if (   boost::math::isnan(calculated)
+            || std::fabs(expected - calculated) > tol) {
+            BOOST_FAIL("failed to reproduce the Lagrange"
+                    " interplation on Chebyshev points"
+                    << "\n    x         : " << x
+                    << "\n    calculated: " << calculated
+                    << "\n    expected  : " << expected);
+        }
+
+        const Real calculatedDeriv = interpl.derivative(x, true);
+        const Real expectedDeriv = std::exp(x)*(std::cos(x) + std::sin(x))
+                / square<Real>()(std::cos(x));
+
+        if (   boost::math::isnan(calculated)
+            || std::fabs(expected - calculated) > tol) {
+            BOOST_FAIL("failed to reproduce the Lagrange"
+                    " interplation derivative on Chebyshev points"
+                    << "\n    x         : " << x
+                    << "\n    calculated: " << calculatedDeriv
+                    << "\n    expected  : " << expectedDeriv);
+        }
+    }
+}
+
+
 test_suite* InterpolationTest::suite() {
     test_suite* suite = BOOST_TEST_SUITE("Interpolation tests");
 
@@ -2013,5 +2185,15 @@ test_suite* InterpolationTest::suite() {
     suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testNoArbSabrInterpolation));
     suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testSabrSingleCases));
     suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testTransformations));
+    suite->add(QUANTLIB_TEST_CASE(
+        &InterpolationTest::testLagrangeInterpolation));
+    suite->add(QUANTLIB_TEST_CASE(
+        &InterpolationTest::testLagrangeInterpolationAtSupportPoint));
+    suite->add(QUANTLIB_TEST_CASE(
+        &InterpolationTest::testLagrangeInterpolationDerivative));
+
+    suite->add(QUANTLIB_TEST_CASE(
+        &InterpolationTest::testLagrangeInterpolationOnChebyshevPoints));
+
     return suite;
 }
diff --git a/test-suite/interpolations.hpp b/test-suite/interpolations.hpp
index d6fca1f1a13..f89cc695057 100644
--- a/test-suite/interpolations.hpp
+++ b/test-suite/interpolations.hpp
@@ -52,6 +52,10 @@ class InterpolationTest {
     static void testNoArbSabrInterpolation();
     static void testSabrSingleCases();
     static void testTransformations();
+    static void testLagrangeInterpolation();
+    static void testLagrangeInterpolationAtSupportPoint();
+    static void testLagrangeInterpolationDerivative();
+    static void testLagrangeInterpolationOnChebyshevPoints();
 
     static boost::unit_test_framework::test_suite* suite();
 };
diff --git a/test-suite/normalclvmodel.cpp b/test-suite/normalclvmodel.cpp
new file mode 100644
index 00000000000..1806c73fde2
--- /dev/null
+++ b/test-suite/normalclvmodel.cpp
@@ -0,0 +1,588 @@
+/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+
+/*
+ Copyright (C) 2016 Klaus Spanderen
+
+ This file is part of QuantLib, a free-software/open-source library
+ for financial quantitative analysts and developers - http://quantlib.org/
+
+ QuantLib is free software: you can redistribute it and/or modify it
+ under the terms of the QuantLib license.  You should have received a
+ copy of the license along with this program; if not, please email
+ <quantlib-dev@lists.sf.net>. The license is also available online at
+ <http://quantlib.org/license.shtml>.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE.  See the license for more details.
+*/
+
+
+#include "utilities.hpp"
+#include "normalclvmodel.hpp"
+
+#include <ql/quotes/simplequote.hpp>
+#include <ql/time/daycounters/actual360.hpp>
+#include <ql/time/daycounters/actualactual.hpp>
+#include <ql/time/daycounters/actual365fixed.hpp>
+#include <ql/time/calendars/nullcalendar.hpp>
+#include <ql/math/statistics/statistics.hpp>
+#include <ql/math/randomnumbers/rngtraits.hpp>
+#include <ql/math/randomnumbers/sobolbrownianbridgersg.hpp>
+#include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp>
+#include <ql/methods/montecarlo/pathgenerator.hpp>
+#include <ql/termstructures/volatility/sabr.hpp>
+#include <ql/termstructures/volatility/equityfx/hestonblackvolsurface.hpp>
+
+#include <ql/experimental/models/normalclvmodel.hpp>
+#include <ql/experimental/finitedifferences/bsmrndcalculator.hpp>
+#include <ql/experimental/finitedifferences/hestonrndcalculator.hpp>
+#include <ql/experimental/finitedifferences/fdornsteinuhlenbeckvanillaengine.hpp>
+#include <ql/experimental/barrieroption/doublebarrieroption.hpp>
+#include <ql/experimental/barrieroption/analyticdoublebarrierbinaryengine.hpp>
+
+#include <boost/bind.hpp>
+#include <boost/make_shared.hpp>
+#include <boost/assign/std/vector.hpp>
+
+using namespace QuantLib;
+using namespace boost::assign;
+using namespace boost::unit_test_framework;
+
+void NormalCLVModelTest::testBSCumlativeDistributionFunction() {
+    BOOST_TEST_MESSAGE("Testing Black-Scholes cumulative distribution function"
+                       " with constant volatility...");
+
+    SavedSettings backup;
+
+    const DayCounter dc = Actual365Fixed();
+    const Date today = Date(22, June, 2016);
+    const Date maturity = today + Period(6, Months);
+
+    const Real s0 = 100;
+    const Real rRate = 0.1;
+    const Real qRate = 0.05;
+    const Volatility vol = 0.25;
+
+    const Handle<Quote> spot(boost::make_shared<SimpleQuote>(s0));
+    const Handle<YieldTermStructure> qTS(flatRate(today, qRate, dc));
+    const Handle<YieldTermStructure> rTS(flatRate(today, rRate, dc));
+    const Handle<BlackVolTermStructure> volTS(flatVol(today, vol, dc));
+
+    const boost::shared_ptr<GeneralizedBlackScholesProcess> bsProcess(
+        boost::make_shared<GeneralizedBlackScholesProcess>(
+            spot, qTS, rTS, volTS));
+    const boost::shared_ptr<OrnsteinUhlenbeckProcess> ouProcess;
+
+    const NormalCLVModel m(
+        bsProcess, ouProcess, std::vector<Date>(), 5);
+    const BSMRNDCalculator rndCalculator(bsProcess);
+
+
+    const Real tol = 1e5*QL_EPSILON;
+    const Time t = dc.yearFraction(today, maturity);
+    for (Real x=10; x < 400; x+=10) {
+        const Real calculated = m.cdf(maturity, x);
+        const Real expected = rndCalculator.cdf(std::log(x), t);
+
+        if (std::fabs(calculated - expected) > tol) {
+            BOOST_FAIL("Failed to reproduce CDF for "
+                       << "\n    strike:     " << x
+                       << "\n    calculated: " << calculated
+                       << "\n    expected:   " << expected);
+        }
+    }
+}
+
+void NormalCLVModelTest::testHestonCumlativeDistributionFunction() {
+    BOOST_TEST_MESSAGE("Testing Heston cumulative distribution function ...");
+
+    SavedSettings backup;
+
+    const DayCounter dc = Actual365Fixed();
+    const Date today = Date(22, June, 2016);
+    const Date maturity = today + Period(1, Years);
+
+    const Real s0 = 100;
+    const Real v0 = 0.01;
+    const Real rRate = 0.1;
+    const Real qRate = 0.05;
+    const Real kappa = 2.0;
+    const Real theta = 0.09;
+    const Real sigma = 0.4;
+    const Real rho = -0.75;
+
+    const Handle<Quote> spot(boost::make_shared<SimpleQuote>(s0));
+    const Handle<YieldTermStructure> qTS(flatRate(today, qRate, dc));
+    const Handle<YieldTermStructure> rTS(flatRate(today, rRate, dc));
+
+    const boost::shared_ptr<HestonProcess> process(
+        boost::make_shared<HestonProcess>(
+            rTS, qTS, spot, v0, kappa, theta, sigma, rho));
+
+    const Handle<BlackVolTermStructure> hestonVolTS(
+        boost::make_shared<HestonBlackVolSurface>(
+            Handle<HestonModel>(boost::make_shared<HestonModel>(process))));
+
+    const NormalCLVModel m(
+        boost::make_shared<GeneralizedBlackScholesProcess>(
+            spot, qTS, rTS, hestonVolTS),
+        boost::shared_ptr<OrnsteinUhlenbeckProcess>(),
+        std::vector<Date>(), 5);
+
+    const HestonRNDCalculator rndCalculator(process);
+
+    const Real tol = 1e-6;
+    const Time t = dc.yearFraction(today, maturity);
+    for (Real x=10; x < 400; x+=25) {
+        const Real calculated = m.cdf(maturity, x);
+        const Real expected = rndCalculator.cdf(std::log(x), t);
+
+        if (std::fabs(calculated - expected) > tol) {
+            BOOST_FAIL("Failed to reproduce CDF for "
+                       << "\n    strike:     " << x
+                       << "\n    calculated: " << calculated
+                       << "\n    expected:   " << expected);
+        }
+    }
+}
+
+
+namespace {
+    class SABRVolTermStructure : public BlackVolatilityTermStructure {
+      public:
+        SABRVolTermStructure(Real alpha, Real beta, Real gamma, Real rho,
+                             Real s0, Real r,
+                             const Date& referenceDate, const DayCounter dc)
+      : BlackVolatilityTermStructure(
+            referenceDate, NullCalendar(), Following, dc),
+        alpha_(alpha), beta_(beta),
+        gamma_(gamma), rho_(rho),
+        s0_(s0), r_(r) { }
+
+        Date maxDate()   const { return Date::maxDate(); }
+        Rate minStrike() const { return 0.0; }
+        Rate maxStrike() const { return QL_MAX_REAL; }
+
+      protected:
+        Volatility blackVolImpl(Time t, Real strike) const {
+            const Real fwd = s0_*std::exp(r_*t);
+            return sabrVolatility(strike, fwd, t, alpha_, beta_, gamma_, rho_);
+        }
+      private:
+        Real alpha_, beta_, gamma_, rho_, s0_, r_;
+    };
+}
+
+void NormalCLVModelTest::testIllustrative1DExample() {
+    BOOST_TEST_MESSAGE("Testing illustrative 1d example of NormalCLVModel ...");
+
+    SavedSettings backup;
+
+    // example taken from:
+    // A. Grzelak, 2015, The CLV Framework -
+    // A Fresh Look at Efficient Pricing with Smile
+    // http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2747541
+
+    const DayCounter dc = Actual360();
+    const Date today = Date(22, June, 2016);
+
+    //SABR
+    const Real beta = 0.5;
+    const Real alpha= 0.2;
+    const Real rho  = -0.9;
+    const Real gamma= 0.2;
+
+    // Ornstein-Uhlenbeck
+    const Real speed = 1.3;
+    const Real level = 0.1;
+    const Real vol   = 0.25;
+    const Real x0    = 1.0;
+
+    const Real s0    = 1.0;
+    const Real rRate = 0.03;
+    const Real qRate = 0.0;
+
+    const Handle<Quote> spot(boost::make_shared<SimpleQuote>(s0));
+    const Handle<YieldTermStructure> qTS(flatRate(today, qRate, dc));
+    const Handle<YieldTermStructure> rTS(flatRate(today, rRate, dc));
+
+    const Handle<BlackVolTermStructure> sabrVol(
+        boost::make_shared<SABRVolTermStructure>(
+            alpha, beta, gamma, rho, s0, rRate, today, dc));
+
+    const boost::shared_ptr<GeneralizedBlackScholesProcess> bsProcess(
+        boost::make_shared<GeneralizedBlackScholesProcess>(
+            spot, qTS, rTS, sabrVol));
+
+    const boost::shared_ptr<OrnsteinUhlenbeckProcess> ouProcess(
+        boost::make_shared<OrnsteinUhlenbeckProcess>(
+            speed, vol, x0, level));
+
+    std::vector<Date> maturityDates;
+    maturityDates += (today + Period(18, Days)),
+        today + Period(90, Days)  , today + Period(180, Days),
+        today + Period(360, Days) , today + Period(720, Days);
+
+    const NormalCLVModel m(bsProcess, ouProcess, maturityDates, 4);
+    const boost::function<Real(Real, Real)> g = m.g();
+
+    // test collocation points in x_ij
+    std::vector<Date> maturities;
+    maturities += maturityDates[0], maturityDates[2], maturityDates[4];
+
+    std::vector<std::vector<Real> > x(3);
+    x[0] += 1.070, 0.984, 0.903, 0.817;
+    x[1] += 0.879, 0.668, 0.472, 0.261;
+    x[2] += 0.528, 0.282, 0.052,-0.194;
+
+    std::vector<std::vector<Real> > s(3);
+    s[0] += 1.104, 1.035, 0.969, 0.895;
+    s[1] += 1.328, 1.122, 0.911, 0.668;
+    s[2] += 1.657, 1.283, 0.854, 0.339;
+
+    std::vector<Real> c;
+    c += 2.3344, 0.7420, -0.7420, -2.3344;
+
+    const Real tol = 0.001;
+    for (Size i=0; i < maturities.size(); ++i) {
+        const Time t = dc.yearFraction(today, maturities[i]);
+
+        for (Size j=0; j < x.front().size(); ++j) {
+            const Real calculatedX = m.collocationPointsX(maturities[i])[j];
+            const Real expectedX = x[i][j];
+
+            if (std::fabs(calculatedX - expectedX) > tol) {
+                BOOST_FAIL("Failed to reproduce collocation x points for "
+                           << "\n    time:       " << maturities[i]
+                           << "\n    j           " << j
+                           << "\n    calculated: " << calculatedX
+                           << "\n    expected:   " << expectedX);
+            }
+
+            const Real calculatedS = m.collocationPointsY(maturities[i])[j];
+            const Real expectedS = s[i][j];
+            if (std::fabs(calculatedS - expectedS) > tol) {
+                BOOST_FAIL("Failed to reproduce collocation s points for "
+                           << "\n    time:       " << maturities[i]
+                           << "\n    j           " << j
+                           << "\n    calculated: " << calculatedS
+                           << "\n    expected:   " << expectedS);
+            }
+
+            const Real expectation
+                = ouProcess->expectation(0.0, ouProcess->x0(), t);
+            const Real stdDeviation
+                = ouProcess->stdDeviation(0.0, ouProcess->x0(), t);
+
+            const Real calculatedG = g(t, expectation + stdDeviation*c[j]);
+            if (std::fabs(calculatedG - expectedS) > tol) {
+                BOOST_FAIL("Failed to reproduce g values "
+                           "at collocation points for "
+                           << "\n    time:       " << maturities[i]
+                           << "\n    j           " << j
+                           << "\n    calculated: " << calculatedG
+                           << "\n    expected:   " << expectedS);
+            }
+        }
+    }
+}
+
+namespace {
+    class NormalCLVModelPayoff : public PlainVanillaPayoff {
+      public:
+        NormalCLVModelPayoff(Option::Type type, Real strike,
+                             const boost::function<Real(Real)> g)
+        : PlainVanillaPayoff(type, strike),
+          g_(g) { }
+
+        Real operator()(Real x) const {
+            return PlainVanillaPayoff::operator()(g_(x));
+        }
+
+      private:
+        const boost::function<Real(Real)> g_;
+    };
+}
+
+void NormalCLVModelTest::testMonteCarloBSOptionPricing() {
+    BOOST_TEST_MESSAGE("Testing Monte-Carlo BS option pricing ...");
+
+    SavedSettings backup;
+
+    const DayCounter dc = Actual365Fixed();
+    const Date today = Date(22, June, 2016);
+    const Date maturity = today + Period(1, Years);
+    const Time t = dc.yearFraction(today, maturity);
+
+    const Real strike = 110;
+    const boost::shared_ptr<StrikedTypePayoff> payoff =
+        boost::make_shared<PlainVanillaPayoff>(Option::Call, strike);
+    const boost::shared_ptr<Exercise> exercise =
+        boost::make_shared<EuropeanExercise>(maturity);
+
+    // Ornstein-Uhlenbeck
+    const Real speed = 2.3;
+    const Real level = 100;
+    const Real sigma = 0.35;
+    const Real x0    = 100.0;
+
+    const Real s0        = x0;
+    const Volatility vol = 0.25;
+    const Real rRate     = 0.10;
+    const Real qRate     = 0.04;
+
+    const Handle<Quote> spot(boost::make_shared<SimpleQuote>(s0));
+    const Handle<YieldTermStructure> qTS(flatRate(today, qRate, dc));
+    const Handle<YieldTermStructure> rTS(flatRate(today, rRate, dc));
+    const Handle<BlackVolTermStructure> vTS(flatVol(today, vol, dc));
+
+    const boost::shared_ptr<GeneralizedBlackScholesProcess> bsProcess(
+        boost::make_shared<GeneralizedBlackScholesProcess>(
+            spot, qTS, rTS, vTS));
+
+    const boost::shared_ptr<OrnsteinUhlenbeckProcess> ouProcess(
+        boost::make_shared<OrnsteinUhlenbeckProcess>(
+            speed, sigma, x0, level));
+
+    std::vector<Date> maturities;
+    maturities += today + Period(6, Months), maturity;
+
+    const NormalCLVModel m(bsProcess, ouProcess, maturities, 8);
+    const boost::function<Real(Real, Real)> g = m.g();
+
+    const Size nSims = 32767;
+    const LowDiscrepancy::rsg_type ld
+        = LowDiscrepancy::make_sequence_generator(1, 23455);
+
+    Statistics stat;
+    for (Size i=0; i < nSims; ++i) {
+        const Real dw = ld.nextSequence().value.front();
+
+        const Real o_t = ouProcess->evolve(0, x0, t, dw);
+        const Real s = g(t, o_t);
+
+        stat.add((*payoff)(s));
+
+    }
+
+    Real calculated = stat.mean() * rTS->discount(maturity);
+
+    VanillaOption option(payoff, exercise);
+    option.setPricingEngine(
+        boost::make_shared<AnalyticEuropeanEngine>(bsProcess));
+    const Real expected = option.NPV();
+
+    const Real tol = 0.01;
+    if (std::fabs(calculated - expected) > tol) {
+        BOOST_FAIL("Failed to reproduce Monte-Carlo vanilla option price "
+                   << "\n    time:       " << maturity
+                   << "\n    strike:     " << strike
+                   << "\n    calculated: " << calculated
+                   << "\n    expected:   " << expected);
+    }
+
+    VanillaOption fdmOption(
+         boost::make_shared<NormalCLVModelPayoff>(
+             payoff->optionType(), payoff->strike(), boost::bind(g, t, _1)),
+         exercise);
+
+    fdmOption.setPricingEngine(
+        boost::make_shared<FdOrnsteinUhlenbeckVanillaEngine>(
+            ouProcess, rTS.currentLink(), 50, 800));
+
+    calculated = fdmOption.NPV();
+    if (std::fabs(calculated - expected) > tol) {
+        BOOST_FAIL("Failed to reproduce FDM vanilla option price "
+                   << "\n    time:       " << maturity
+                   << "\n    strike:     " << strike
+                   << "\n    calculated: " << calculated
+                   << "\n    expected:   " << expected);
+    }
+}
+
+void NormalCLVModelTest::testMoustacheGraph() {
+    BOOST_TEST_MESSAGE(
+        "Testing Double no Touch pricing with Normal CLV Model...");
+
+    SavedSettings backup;
+
+    /*
+     The comparison of Black-Scholes and normal CLV prices is derived
+     from figure 8.8 in Iain J. Clark's book,
+     Foreign Exchange Option Pricing: A Practitioner’s Guide
+    */
+
+    const DayCounter dc = ActualActual();
+    const Date todaysDate(5, Aug, 2016);
+    const Date maturityDate = todaysDate + Period(1, Years);
+    const Time maturityTime = dc.yearFraction(todaysDate, maturityDate);
+
+    Settings::instance().evaluationDate() = todaysDate;
+
+    const Real s0 = 100;
+    const Handle<Quote> spot(boost::make_shared<SimpleQuote>(s0));
+    const Rate r = 0.02;
+    const Rate q = 0.01;
+
+    // parameter of the "calibrated" Heston model
+    const Real kappa =   1.0;
+    const Real theta =   0.06;
+    const Real rho   =  -0.8;
+    const Real sigma =   0.8;
+    const Real v0    =   0.09;
+
+    const Handle<YieldTermStructure> rTS(flatRate(r, dc));
+    const Handle<YieldTermStructure> qTS(flatRate(q, dc));
+
+    const boost::shared_ptr<HestonModel> hestonModel(
+        boost::make_shared<HestonModel>(
+            boost::make_shared<HestonProcess>(
+                rTS, qTS, spot, v0, kappa, theta, sigma, rho)));
+
+    const Handle<BlackVolTermStructure> vTS(
+        boost::make_shared<HestonBlackVolSurface>(
+            Handle<HestonModel>(hestonModel)));
+
+    const boost::shared_ptr<GeneralizedBlackScholesProcess> bsProcess =
+        boost::make_shared<GeneralizedBlackScholesProcess>(
+            spot, qTS, rTS, vTS);
+
+    // Ornstein-Uhlenbeck
+    const Real speed   = -0.80;
+    const Real level   = 100;
+    const Real sigmaOU = 0.15;
+    const Real x0      = 100;
+
+    const boost::shared_ptr<OrnsteinUhlenbeckProcess> ouProcess(
+        boost::make_shared<OrnsteinUhlenbeckProcess>(
+            speed, sigmaOU, x0, level));
+
+    const boost::shared_ptr<Exercise> europeanExercise(
+        boost::make_shared<EuropeanExercise>(maturityDate));
+
+    VanillaOption vanillaOption(
+        boost::make_shared<PlainVanillaPayoff>(Option::Call, s0),
+        europeanExercise);
+
+    vanillaOption.setPricingEngine(
+        boost::make_shared<AnalyticHestonEngine>(hestonModel));
+
+    const Volatility atmVol = vanillaOption.impliedVolatility(
+        vanillaOption.NPV(),
+        boost::make_shared<GeneralizedBlackScholesProcess>(spot, qTS, rTS,
+            Handle<BlackVolTermStructure>(flatVol(std::sqrt(theta), dc))));
+
+    const boost::shared_ptr<PricingEngine> analyticEngine(
+        boost::make_shared<AnalyticDoubleBarrierBinaryEngine>(
+            boost::make_shared<GeneralizedBlackScholesProcess>(
+                spot, qTS, rTS,
+                Handle<BlackVolTermStructure>(flatVol(atmVol, dc)))));
+
+
+    std::vector<Date> maturities(1, todaysDate + Period(2, Weeks));
+    while (maturities.back() < maturityDate)
+        maturities.push_back(maturities.back() + Period(2, Weeks));
+
+    const NormalCLVModel m(bsProcess, ouProcess, maturities, 8);
+    const boost::function<Real(Real, Real)> g = m.g();
+
+    const Size n = 18;
+    Array barrier_lo(n), barrier_hi(n), bsNPV(n);
+
+    const boost::shared_ptr<CashOrNothingPayoff> payoff =
+        boost::make_shared<CashOrNothingPayoff>(Option::Call, 0.0, 1.0);
+
+    for (Size i=0; i < n; ++i) {
+        const Real dist = 10.0+5.0*i;
+
+        barrier_lo[i] = std::max(s0 - dist, 1e-2);
+        barrier_hi[i] = s0 + dist;
+        DoubleBarrierOption doubleBarrier(
+            DoubleBarrier::KnockOut, barrier_lo[i], barrier_hi[i], 0.0,
+            payoff,
+            europeanExercise);
+
+        doubleBarrier.setPricingEngine(analyticEngine);
+        bsNPV[i] = doubleBarrier.NPV();
+    }
+
+    typedef SobolBrownianBridgeRsg rsg_type;
+    typedef PathGenerator<rsg_type>::sample_type sample_type;
+
+    const Size factors = 1;
+    const Size tSteps = 200;
+    const TimeGrid grid(maturityTime, tSteps);
+
+    const boost::shared_ptr<PathGenerator<rsg_type> > pathGenerator =
+        boost::make_shared<PathGenerator<rsg_type> >(
+            ouProcess, grid, rsg_type(factors, tSteps), false);
+
+    const Size nSims = 400000;
+    std::vector<GeneralStatistics> stats(n);
+    const DiscountFactor df = rTS->discount(maturityDate);
+
+    for (Size i=0; i < nSims; ++i) {
+        std::vector<bool> touch(n, false);
+
+        const sample_type& path = pathGenerator->next();
+
+        Real s;
+        for (Size j=1; j <= tSteps; ++j) {
+            const Time t = grid.at(j);
+            s = g(t, path.value.at(j));
+
+            for (Size u=0; u < n; ++u) {
+                if (s <= barrier_lo[u] || s >= barrier_hi[u]) {
+                    touch[u] = true;
+                }
+            }
+        }
+        for (Size u=0; u < n; ++u) {
+            if (touch[u]) {
+                stats[u].add(0.0);
+            }
+            else {
+                stats[u].add(df*payoff->operator()(s));
+            }
+        }
+    }
+
+    const Real expected[] = {
+        0.00907689,  0.0892242,  0.139655,   0.111603,   0.059872,
+        0.016652,   -0.0094134, -0.0207528, -0.0230873, -0.0201951,
+        -0.0160556, -0.0116482, -0.00741535,-0.00355287,-0.000165914,
+        0.00493464,  0.00342315, 0.0023405 };
+
+    const Real tol = 1e-5;
+    for (Size u=0; u < n; ++u) {
+        const Real calculated = stats[u].mean() - bsNPV[u];
+
+        if (std::fabs(calculated - expected[u]) > tol) {
+            BOOST_FAIL("Failed to reproduce Double no Touch prices"
+                   << "\n    time:          " << maturityDate
+                   << "\n    barrier lower: " << barrier_lo[u]
+                   << "\n    barrier high:  " << barrier_hi[u]
+                   << "\n    calculated:    " << calculated
+                   << "\n    expected:      " << expected[u]);
+        }
+    }
+}
+
+
+test_suite* NormalCLVModelTest::experimental() {
+    test_suite* suite = BOOST_TEST_SUITE("NormalCLVModel tests");
+
+    suite->add(QUANTLIB_TEST_CASE(
+        &NormalCLVModelTest::testBSCumlativeDistributionFunction));
+    suite->add(QUANTLIB_TEST_CASE(
+        &NormalCLVModelTest::testHestonCumlativeDistributionFunction));
+    suite->add(QUANTLIB_TEST_CASE(
+        &NormalCLVModelTest::testIllustrative1DExample));
+    suite->add(QUANTLIB_TEST_CASE(
+            NormalCLVModelTest::testMonteCarloBSOptionPricing));
+    suite->add(QUANTLIB_TEST_CASE(
+            NormalCLVModelTest::testMoustacheGraph));
+
+    return suite;
+}
diff --git a/test-suite/normalclvmodel.hpp b/test-suite/normalclvmodel.hpp
new file mode 100644
index 00000000000..069b7fa2e97
--- /dev/null
+++ b/test-suite/normalclvmodel.hpp
@@ -0,0 +1,38 @@
+/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+
+/*
+ Copyright (C) 2016 Klaus Spanderen
+
+ This file is part of QuantLib, a free-software/open-source library
+ for financial quantitative analysts and developers - http://quantlib.org/
+
+ QuantLib is free software: you can redistribute it and/or modify it
+ under the terms of the QuantLib license.  You should have received a
+ copy of the license along with this program; if not, please email
+ <quantlib-dev@lists.sf.net>. The license is also available online at
+ <http://quantlib.org/license.shtml>.
+
+ This program is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ FOR A PARTICULAR PURPOSE.  See the license for more details.
+*/
+
+#ifndef quantlib_test_normal_clv_model_hpp
+#define quantlib_test_normal_clv_model_hpp
+
+#include <boost/test/unit_test.hpp>
+
+
+class NormalCLVModelTest {
+   public:
+    static void testBSCumlativeDistributionFunction();
+    static void testHestonCumlativeDistributionFunction();
+    static void testIllustrative1DExample();
+    static void testMonteCarloBSOptionPricing();
+    static void testMoustacheGraph();
+
+    static boost::unit_test_framework::test_suite* experimental();
+};
+
+
+#endif
diff --git a/test-suite/quantlibtestsuite.cpp b/test-suite/quantlibtestsuite.cpp
index 25d3a99f902..10b326d0cd1 100644
--- a/test-suite/quantlibtestsuite.cpp
+++ b/test-suite/quantlibtestsuite.cpp
@@ -152,6 +152,7 @@
 #include "mersennetwister.hpp"
 #include "money.hpp"
 #include "noarbsabr.hpp"
+#include "normalclvmodel.hpp"
 #include "nthtodefault.hpp"
 #include "numericaldifferentiation.hpp"
 #include "observable.hpp"
@@ -457,6 +458,7 @@ test_suite* init_unit_test_suite(int, char* []) {
     test->add(InflationVolTest::suite());
     test->add(MargrabeOptionTest::suite());
     test->add(NoArbSabrTest::suite());
+    test->add(NormalCLVModelTest::experimental());
     test->add(NthToDefaultTest::suite());
     test->add(NumericalDifferentiationTest::suite());
     test->add(PagodaOptionTest::suite());
diff --git a/test-suite/testsuite.vcxproj b/test-suite/testsuite.vcxproj
index 6b0ef3ba887..b97b1a5eadb 100644
--- a/test-suite/testsuite.vcxproj
+++ b/test-suite/testsuite.vcxproj
@@ -719,6 +719,7 @@
     <ClCompile Include="mersennetwister.cpp" />
     <ClCompile Include="money.cpp" />
     <ClCompile Include="noarbsabr.cpp" />
+    <ClCompile Include="normalclvmodel.cpp" />
     <ClCompile Include="nthtodefault.cpp" />
     <ClCompile Include="numericaldifferentiation.cpp" />
     <ClCompile Include="observable.cpp" />
@@ -860,6 +861,7 @@
     <ClInclude Include="mersennetwister.hpp" />
     <ClInclude Include="money.hpp" />
     <ClInclude Include="noarbsabr.hpp" />
+    <ClInclude Include="normalclvmodel.hpp" />
     <ClInclude Include="nthtodefault.hpp" />
     <ClInclude Include="numericaldifferentiation.hpp" />
     <ClInclude Include="observable.hpp" />
@@ -871,7 +873,7 @@
     <ClInclude Include="pagodaoption.hpp" />
     <ClInclude Include="partialtimebarrieroption.hpp" />
     <ClInclude Include="pathgenerator.hpp" />
-    <ClInclude Include="paralleltestrunner.hpp" />    
+    <ClInclude Include="paralleltestrunner.hpp" />
     <ClInclude Include="period.hpp" />
     <ClInclude Include="piecewiseyieldcurve.hpp" />
     <ClInclude Include="piecewisezerospreadedtermstructure.hpp" />
diff --git a/test-suite/testsuite.vcxproj.filters b/test-suite/testsuite.vcxproj.filters
index 0c44c72f643..f41ac79961a 100644
--- a/test-suite/testsuite.vcxproj.filters
+++ b/test-suite/testsuite.vcxproj.filters
@@ -426,6 +426,9 @@
     <ClCompile Include="observable.cpp">
       <Filter>Source Files</Filter>
     </ClCompile>
+    <ClCompile Include="normalclvmodel.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="americanoption.hpp">
@@ -848,5 +851,8 @@
     <ClInclude Include="observable.hpp">
       <Filter>Header Files</Filter>
     </ClInclude>
+    <ClInclude Include="normalclvmodel.hpp">
+      <Filter>Header Files</Filter>
+    </ClInclude>
   </ItemGroup>
 </Project>