Рада всех приветствовать в этот субботний день. Вернулась к торговле, но как-то скучно и пустовато стало в опционных стаканах. И вот — так совпало © — встретился мне прекрасный вдумчивый автор со статьями на тему дельта-хеджирования, и решила я перевести его статью и зашить в блокнотик Jupyter, чем и делюсь с сообществом. Статья состоит из пяти частей. Первые две части я в силу очевидности изложенного в них опустила. Начала перевод с части 3. Так что будут еще два блокнота по части 4 (про ошибку репликации и частоту хеджа) и 5 (про последствия хеджирования по «неправильной» дельте). Не переключайте канал

Марк Джеймисон «Как дельтахеджировать опционы» Блокнот № 1

Может быть лучик просвещения даст сомневающимся последний ключик для начала торговли!

Чтобы утащить блокнот себе, жмите Файл — Сохранить копию на Диске — и сохраняйте у себя на диске. Сохраненная копия будет редактируемой.

Всем добра и профитов

В комментарии к посту пользователь Мальчик buybuy попросил перевернуть те графики, которые были на полшестого, а также прислать полные данные по эксперименту.
Никаких проблем. Однако, чтобы зрителям не было скучно, в текущем посте рассмотрим картинки не одного, а нескольких индикаторов:

1) на основе СЛАУ с одной переменной

indicator = d1 * (d1 / d2)

2) на основе СЛАУ с двумя переменными (см. предыдущий пост)

indicator = d1 * (d1 * d4 - d2 * d3) + d2 * (d2 * d2 - d1 * d3)

3) на основе СЛАУ с тремя переменными.

indicator = d1 * (-d1 * d4 * d6 + d1 * d5 * d5 + d2 * d3 * d6 - d2 * d4 * d5 + d3 * d3 * (-d5) + d3 * d4 * d4) /<br />  (-d2 * d4 * d6 + d2 * d5 * d5 + d3 * d3 * d6 - 2 * d3 * d4 * d5 + d4 * d4 * d4) +<br /><br />  d2 * (-d1 * d3 * d6 + d1 * d4 * d5 + d2 * d2 * d6 - d2 * d3 * d5 + d4 * d4 * d4) /<br />  (d2 * d4 * d6 - d2 * d5 * d5 - d3 * d3 * d6 + 2 * d3 * d4 * d5 - d4 * d4 * d4) +<br /><br />  d3 * (-d1 * d3 * d5 + d1 * d4 * d4 + d2 * d2 * d5 - 2 * d2 * d3 * d4 + d3 * d3 * d3) /<br />  (-d2 * d4 * d6 + d2 * d5 * d5 + d3 * d3 * d6 - 2 * d3 * d4 * d5 + d4 * d4 * d4)

4) на основе СЛАУ с 4 переменными (формулу смотреть в исходниках).
Пока, думаю, этого достаточно (если, конечно, Русский ВПК не укажет дальнейшее направление).

( Читать дальше )

Обращение ко всем потенциальным клиентам АО «ФИНАМ» с целью уберечь их от возможных потерь!!!

Трагедия и полное разочарование «профессионализмом» сотрудников АО «ФИНАМ» в Екатеринбурге!

В апреле 2017 года я хотел открыть вклад под 8% годовых в рублях и положить денежные средства в АО «Банк ФИНАМ» в Екатеринбурге, однако был введен в заблуждение сотрудниками АО «ФИНАМ», которые гарантировали мне гораздо больший процентный доход на ИИС в АО «ФИНАМ», чем на вкладе в АО «Банк ФИНАМ», в частности речь шла об одной из стратегий с «гарантированной» доходностью до 20% годовых в рублях (ОФЗ).

Спустя 2 года сотрудники АО «ФИНАМ вынудили меня продать ОФЗ под предлогом их якобы падающей доходности и перевести мои денежные средства на «самую эффективную», с их слов, стратегию с «гарантированной» доходностью до 20% годовых в рублях.

За 3 года доходность моих денежных средств на ИИС в АО „ФИНАМ“ провалилась до –40%, но и это еще не предел. Далее сотрудники АО „ФИНАМ“ предложили мне перевести средства на «самую консервативную» валютную стратегию, но, несмотря на это мои денежные средства продолжили уменьшаться.

( Читать дальше )

Вчера написал пост — сколько нужно учиться трейдингу? Типа, спросил, т.к. ранее многие утверждали, что учиться нужно 5, а то и все 10 лет.
Ну, и попросил рассказать чему же конкретно в трейдинге, а заодно и в инвестициях нужно так долго учиться? Ну, это же просто, перечисли дисциплины, типа, физика, математика, химия, история, биология. Как в школе. Ну, или хотя бы разделы — нужна, скажем, таблица умножения и тригонометрия, и без этого ну никак.
В итоге, нужно чему-то учиться лет этак 5, лучше 10, чему и как никому неизвестно, но учиться, учиться,… Тайные знания, видимо.)
В общем, как я и предполагал нет никаких тайных знаний, обучаться, в общем, нечему, да и сам предмет обучения просто отсутствует. И вообще, в трейдинге и инвестициях вообще все просто, понятно и доступно. Не получается? — ну, ребята, не все в состоянии даже школу нормально закончить, а уж в институте к 6-му курсу оставалось полгруппы. Стало быть, не ваше, и только.
С трейдингом вообще все просто, и долго на нем останавливаться не буду.

( Читать дальше )

Смотрю люди продолжают нести деньги мошенникам в надежде урвать грааль.


Истина же заключается в том, что рабочую систему, позволяющую масштабировать прибыль вам ни кто ни когда не продаст, ни за какие деньги.

Я вот тоже долго искал комфортную ТС на срочном рынке.

В итоге остановился на шорте фьюча лесенкой под наличный бакс (евро).

Что имеем.

1. 10 лотов наличного евро, купленных по 70 руб в 20 году.

2. 10 лотов фьючерса евро/рубля продаваемых лесенкой со спредом 2 руб.

3. Фикс прибыли от 2000 руб.

4. Ежедневно выплачиваемое контанго по фьючерсу, в среднем 1200-1500 пунктов каждые 3 месяца.

5. ГО 10 тыс руб и ОФЗ.

Риски

1. Уход курса выше 95 руб — тупо закрою срочку наличным евро.

2. Проблемы с брокером, отсутствие времени - тупо закрою срочку наличным евро.


Итоги за год -
по наличному евро — 160 тыс руб прибыли
по срочке 170 тыс прибыли.

Общая доходность, с учетом размещенных ден средств (10 лотов евро) 750 тыс руб, около 45%.

( Читать дальше )

Как им пользоваться:
в папке QUIK\LuaIndicators создаем текстовый файл с содержимым кода индикатора
добавляем индикатор на график
задаем параметр чувствиельности инлдикатора — Процент
индикатор показывает уровни, в которых находится тренд
если тренд еще не сформировался, то границы канала могут расширяться по мере движения цены
если в канале уже есть локальные вершины, то при следующем походе индикатора к границе канал возмен разворот тренда или коррекции в рамках данного канала

--[[
индикатор торговой системы QUIK: 
канал на последнем зигзаге
индикатор строит канал по зигзагу, показывает слом тренда при выходе цены за границы канала
параметры:
  Procent - величина зигзага в %
  k - коэффизиент ширины канала
  Ext - =0-канал стоится по средним =1-канал стоится по экстремуму
Как им пользоваться:
в папке QUIK\LuaIndicators создаем текстовый файл с содержимым кода индикатора
добавляем индикатор на график
задаем параметр чувствиельности инлдикатора — Процент
индикатор показывает уровни, в которых находится тренд
если тренд еще не сформировался, то границы канала могут расширяться по мере движения цены
если в канале уже есть локальные вершины, то при следующем походе индикатора к границе канал возмен разворот тренда или коррекции в рамках данного канала  
--]]
Settings=              
        {                          
            Name = "Zigzag_channel_fast",   -- название индикатора
            Procent=5.0,               -- величина зигзага в % 
            k=2,			         -- коэффизиент ширины канала
            Ext=1, -- =0-канал стоится по средним =1-канал стоится по экстремуму
            shif=0,			
            line=                                     
                {                               
                    {  
                        Name = "zigzagline",
                        Type =TYPE_LINE,
                        Width = 1,
                        Color = RGB(120,90, 140)
                    },
					{  
                        Name = "cur1",
                        Type =TYPE_LINE,
                        Width = 1,
                        Color = RGB(0,0, 0)
                    },
					{  
                        Name = "cur2",
                        Type =TYPE_LINE,
                        Width = 2,
                        Color = RGB(0,200, 200)
                    },
					{  
                        Name = "cur3",
                        Type =TYPE_LINE,
                        Width = 2,
                        Color = RGB(200,200, 0)
                    },
					{  
                        Name = "cur4",
                        Type =TYPE_LINE,
                        Width = 1,
                        Color = RGB(0,0, 0)
                    },
					{  
                        Name = "cur5",
                        Type =TYPE_LINE,
                        Width = 2,
                        Color = RGB(0,0, 255)
                    },
					{  
                        Name = "cur6",
                        Type =TYPE_LINE,
                        Width = 2,
                        Color = RGB(255,0, 0)
                    },
					{  
                        Name = "cur7",
                        Type =TYPE_LINE,
                        Width = 2,
                        Color = RGB(255,0, 0)
                    },
					{  
                        Name = "cur8",
                        Type =TYPE_LINE,
                        Width = 2,
                        Color = RGB(255,0, 0)
                    },
					{  
                        Name = "cur9",
                        Type =TYPE_LINE,
                        Width = 2,
                        Color = RGB(255,0, 0)
                    },
					{  
                        Name = "cur10",
                        Type =TYPE_LINE,
                        Width = 2,
                        Color = RGB(255,0, 0)
                    },
					{  
                        Name = "cur11",
                        Type =TYPE_LINE,
                        Width = 2,
                        Color = RGB(255,0, 0)
                    }
                }
       }


function Init()

  y1 = nil
  y2 = nil
  y3 = nil  
  x1 = 1
  x2 = 1
  x3 = 1
  
  oldindex = 1
  
  a_01 = 0
  b_01 = 0
  a_02 = 0
  b_02 = 0
  
  vminus_01 = 0
  vplus_01 = 0  
  vminus_02 = 0
  vplus_02 = 0  
  vFrom_01 = 1
  vTo_01 = 1
  vFrom_02 = 1
  vTo_02 = 1
  
  vminus = 0
  vplus = 0
  vminus_n = 0
  vplus_n = 0  
  
  ret = 0
  
  last_ind = 1

  v1=nil
  v2=nil
  v3=nil
  v4=nil
  v5=nil
  v6=nil
  v7=nil
  v8=nil
  v9=nil
  v10=nil
  
  return 10
end

function recalcline(num, index, rec)

 if index == Size()-sh and rec == 0 then
   for i = 1, index do
     if m == 1 then 		
	   SetValue(i, 3, nil)
	   SetValue(i, 4, nil)
	 end
     if m == 2 then 		
	   SetValue(i, 6, nil)
	   SetValue(i, 7, nil)
	 end 	 
   end 
 end 

 if index == Size()-sh or rec == 1 then
	
  if num == 1 then 
	  vFrom = 0 --[[ --]] 
	  vTo = 0
	if x3 < x2 then 
	    vFrom = x3
	    vTo = x2
	  else
	    vFrom = x2
	    vTo = x3
	  end 
  end
  if num == 2 then 
	  vFrom = 0
	  vTo = 0
	  if x1 < x2 then 
	    vFrom = x1
	    vTo = x2
	  else
	    vFrom = x2
	    vTo = x1
	  end 
  end    
  
      if vFrom < 1 then
	    vFrom = 1
	  end 	
  
	  a1 = 0
	  a2 = 0
	  a3 = 0
	  a4 = 0
	  n = 0
	  a = 0
	  b = 0
      for i = vFrom, vTo do
        n=n+1	  
	    a1 = a1+i*C(i)
	    a2 = a2+i
	    a3 = a3+C(i)
	    a4 = a4+i*i	
	  end
	  if((n*a4 - a2*a2) ~= 0) then	  
	  
	  -- if oldindex ~= Size() then
	   a = (n*a1 - a2*a3)/(n*a4 - a2*a2)
	   b = (a3 - a*a2)/n

	   vminus = 0
	   vplus = 0
	   vminus_n = 0
	   vplus_n = 0
       for i = vFrom, index do
        y = a*i + b		
		if i>0  then 
		if C(i) < y and i <= vTo then
		  if ex == 0 then 
	        vminus = vminus + y - C(i)
            vminus_n = vminus_n + 1	
		  else 
		    if i == vFrom then 
			  vminus = 0
			end
            if y-C(i) > vminus then
              vminus = y-C(i) 			
            end				
		  end 
        end	
        end		
		if i>0  then 
		if C(i) > y and i <= vTo then
		  if ex == 0 then 
	        vplus = vplus + C(i) - y
	        vplus_n = vplus_n + 1
		  else 
		    if i == vFrom then 
			  vplus = 0
			end		  
            if C(i)-y > vplus then
              vplus = C(i)-y 			
            end			
		  end 
		end 
		end  
        m = 0
        if C(vTo) >	C(vFrom) and y < C(vTo)+0.1*(C(vTo)-C(vFrom)) and y > C(vFrom)-0.1*(C(vTo)-C(vFrom)) then
          m = 1
		end 	
        if C(vTo) <	C(vFrom) and y < C(vFrom)+0.1*(C(vFrom)-C(vTo)) and y > C(vTo)-0.1*(C(vFrom)-C(vTo)) then
          m = 1
		end 
        if m == 1 then 		
	      --SetValue(i, 2, y)
        end 		
	   end
	   --end
	  
	   if ex == 0 then  
	    if vminus_n ~= 0 then
	     vminus = kk*vminus/vminus_n          
	    end 
	    if vplus_n ~= 0 then
	     vplus = kk*vplus/vplus_n          
	    end 
	   end 	  
      end

  if num == 1 then 
    a_01 = a
    b_01 = b 
    vminus_01 = vminus
    vplus_01 = vplus
	vFrom_01 = vFrom
    vTo_01 = vTo
  end
  if num == 2 then 
    a_02 = a
    b_02 = b  
    vminus_02 = vminus
    vplus_02 = vplus
	vFrom_02 = vFrom
    vTo_02 = vTo	
  end  
  return 1 
 end
 
 return 0
end 

function levels(x, y, num, index)

  if y ~= nil and x ~= 1 then 
  
    for i = 1, index do	
	  if i > x then 
	    SetValue(i, num, y)
      else 		
	    SetValue(i, num, nil)
	  end 
    end
	
  end 

end


function OnCalculate(index)
  
  de = Settings.Procent
  kk = Settings.k 
  ex = Settings.Ext
  sh = Settings.shif
  
  v1=nil
  v2=nil
  v3=nil
  v4=nil
  v5=nil
  v6=nil
  v7=nil  
  v8=nil
  v9=nil
  v10=nil  
  
   
  if index > Size()-sh then  
    --return nil, nil, nil, nil, nil, nil, nil 	 
    return v1, v2, v3, v4, v5, v6, v7, v8, v9, v10 
  end 
  
  vl = C(index)
  if vl~=nil then 
   if index == 1 then 
	y1 = vl
    y2 = vl
   else   
	  if C(index) > y1*(1+de/100) and y1 < y2 then 
	    x3 = x2
	    x2 = x1
		y3 = y2
	    y2 = y1	
	    x1 = index 
	    y1 = C(index)	
        ret = recalcline(1, index, 0)
        ret = recalcline(2, index, 0)   
        levels(x1, y1, 8, index)		
		levels(x2, y2, 9, index)
		levels(x3, y3, 10, index)
	  end 	
	  if C(index) > y1 and C(index) > y2 then 
	    x1 = index 
	    y1 = C(index)	  			  
	  end 
	  	  		
	  if C(index) < y1*(1-de/100) and y1 > y2 then 
	    x3 = x2
	    x2 = x1
		y3 = y2		
	    y2 = y1
	    x1 = index 
	    y1 = C(index)	  	
        ret = recalcline(1, index, 0)
        ret = recalcline(2, index, 0)
        levels(x1, y1, 8, index)		
		levels(x2, y2, 9, index)
		levels(x3, y3, 10, index)
	  end 	
	  if C(index) < y1 and C(index) < y2 then 
	    x1 = index 
	    y1 = C(index)	  			  
	  end 	  	  		
	end
  end 	
	
    --[[ 
  if x1 ~= index then 
    curfrom = x1
	curto = index
  else 
    curfrom = x2
	curto = x1
  end 
 
  if curto ~= curfrom and curfrom ~= nil and curto ~= nil then 
    if C(curto) ~= nil and C(curfrom) ~= nil then 
      k = (C(curto)- C(curfrom))/(curto- curfrom)  
      for i = curfrom, index  do
        curv = i*k + C(curto) - curto*k  		          
	    SetValue(i, 1, curv)
      end   	
	end 
  end 
  

 	--]]

	if index == Size()-sh and ret == 0 then
        lpr = C(Size()-sh)
        ret = recalcline(1, index, 1)
        ret = recalcline(2, index, 1)	
		
		if y1 ~= nil then 		
		 if y1 >	lpr*0.9 and y1 < lpr*1.1 then 
		  levels(x1, y1, 8, index)
         end 
        end 		
        if y2 ~= nil then 				
		 if y2 > lpr*0.9 and y2 < lpr*1.1 then 		
		  levels(x2, y2, 9, index)
         end 
		end 
        if y3 ~= nil then 		
		 if y3 > lpr*0.9 and y3 < lpr*1.1 then 		
		  levels(x3, y3, 10, index)
         end 		  
		end 
    
    end 	
	
	
	
	if index == Size()-sh  then  --- and index ~= last_ind
		  
	   if vTo_01 - vFrom_01 > 500 then
         vFrom_01 = vTo_01-500 
       end 	  
	   
       for i = 1, index do		   		
	     SetValue(i, 3, nil)	
         SetValue(i, 4, nil)	 
	     SetValue(i, 6, nil)	
         SetValue(i, 7, nil)	
		-- SetValue(i, 8, nil)	 
	   --  SetValue(i, 9, nil)	
        -- SetValue(i, 10, nil)	
       end	   
		  
       for i = vFrom_01, index do
        y = a_01*i + b_01 + vplus_01		
        m = 0
        if C(vTo_01) >	C(vFrom_01) and y < C(vTo_01)+0.3*(C(vTo_01)-C(vFrom_01)) and y > C(vFrom_01)-0.3*(C(vTo_01)-C(vFrom_01)) then
          m = 1
		end 	
        if C(vTo_01) <	C(vFrom_01) and y < C(vFrom_01)+0.3*(C(vFrom_01)-C(vTo_01)) and y > C(vTo_01)-0.3*(C(vFrom_01)-C(vTo_01)) then
          m = 1
		end 
        if m == 1 then 		
	      SetValue(i, 3, y)
		  if i == index then 
		    v3 = y 
		  end 
        end 	
		
        y = a_01*i + b_01 - vminus_01		
        m = 0
        if C(vTo_01) >	C(vFrom_01) and y < C(vTo_01)+0.3*(C(vTo_01)-C(vFrom_01)) and y > C(vFrom_01)-0.3*(C(vTo_01)-C(vFrom_01)) then
          m = 1
		end 	
        if C(vTo_01) <	C(vFrom_01) and y < C(vFrom_01)+0.3*(C(vFrom_01)-C(vTo_01)) and y > C(vTo_01)-0.3*(C(vFrom_01)-C(vTo_01)) then
          m = 1
		end 
        if m == 1 then 		
	      SetValue(i, 4, y)
		  if i == index then 
		    v4 = y 
		  end 
        end 			
	   end	  
	    
	   if vTo_02 - vFrom_02 > 500 then
         vFrom_02 = vTo_02-500  
       end 	 
	  
       for i = vFrom_02, index do
        y = a_02*i + b_02 + vplus_02		
        m = 0
        if C(vTo_02) >	C(vFrom_02) and y < C(vTo_02)+0.1*(C(vTo_02)-C(vFrom_02)) and y > C(vFrom_02)-0.1*(C(vTo_02)-C(vFrom_02)) then
          m = 1
		end 	
        if C(vTo_02) <	C(vFrom_02) and y < C(vFrom_02)+0.1*(C(vFrom_02)-C(vTo_02)) and y > C(vTo_02)-0.1*(C(vFrom_02)-C(vTo_02)) then
          m = 1
		end 
        if m == 1 then 		
	      SetValue(i, 6, y)
		  if i == index then 
		    v6 = y 
		  end 
        end 	
		
        y = a_02*i + b_02 - vminus_02		
        m = 0
        if C(vTo_02) >	C(vFrom_02) and y < C(vTo_02)+0.1*(C(vTo_02)-C(vFrom_02)) and y > C(vFrom_02)-0.1*(C(vTo_02)-C(vFrom_02)) then
          m = 1
		end 	
        if C(vTo_02) <	C(vFrom_02) and y < C(vFrom_02)+0.1*(C(vFrom_02)-C(vTo_02)) and y > C(vTo_02)-0.1*(C(vFrom_02)-C(vTo_02)) then
          m = 1
		end 
        if m == 1 then 		
	      SetValue(i, 7, y)
		  if i == index then 
		    v7 = y 
		  end 
        end 			
	   end	
	  
	  last_ind = index 
	end 
  
 
    return v1, v2, v3, v4, v5, v6, v7, v8, v9, v10 

end

--[[
параметры: 
Procent - процент зигзага 
--]]
Settings={
Name="ZIGZAG_Templ",
Procent=2,
    line=                                     
                {  
					{  
                        Name = "cur1",
                        Type =TYPE_LINE,
                        Width = 2,
                        Color = RGB(0,0, 0)
                    }				
                }
}

function Init()
  
  y1 = nil
  y2 = nil
  x1 = 1
  x2 = 1
       
  return 1
  
end

function OnCalculate(index)

  de = Settings.Procent

  vl = C(index)
  if index == 1 then 
	y1 = vl
    y2 = vl
  else   
	  if C(index) > y1*(1+de/100) and y1 < y2 then 
	    x2 = x1
	    y2 = y1	
	    x1 = index 
	    y1 = C(index)		
      else 
	    if C(index) > y1 and y1 >= y2 
		then 
	      x1 = index 
	      y1 = C(index)	  			  
	    end 		
	  end 	

	  	  		
	  if C(index) < y1*(1-de/100) and y1 > y2 then 
	    x2 = x1
	    y2 = y1
	    x1 = index 
	    y1 = C(index)				
      else 
	    if C(index) < y1 and y1 <= y2 
		then 
	      x1 = index 
	      y1 = C(index)	  			  
	    end 		
	  end 	
	  	  		
	end 	
  
  if x1 ~= index then 
    curfrom = x1
	curto = index
  else 
    curfrom = x2
	curto = x1
  end 

  if curto ~= curfrom and curfrom ~= nil and curto ~= nil then 
    if C(curto) ~= nil and C(curfrom) ~= nil then 
      k = (C(curto)- C(curfrom))/(curto- curfrom)  
      for i = curfrom, index  do
        curv = i*k + C(curto) - curto*k  		          
	    SetValue(i, 1, curv)
      end   	
	end 
  end 
  
  return vl
 
  
end

После введения утренней торговой сессии проблема автоматического запуска торгового ПО стала особенно актуальна.
Хорошее решение предложил Евгений Логунов  в своей статье «Простой автологин за 5 минут».  Мы предложим аналогичное решение для КВИК на С++.

Задача очень простая — в 7:00 пробудить ПК с помощью планировщика заданий Windows, запустить несколько терминалов QUIK, и в каждом из них запустить торговых роботов, чтобы полностью освободить владельца всего этого счастья от физических и психических нагрузок, плохо влияющих на питание и здоровый образ жизни.


Итак, первое, что нам необходимо будет сделать это Автологин. Штука достаточно простая, учитывая то, что после запуска терминала он автоматически выдает окно приветствия. Нам нужно только дождаться появления этого окна, получить первое вводное поле (логин), второе вводное поле (пароль), вбить туда нужные значения и нажать на первую дочернюю кнопку этого окна: «Вход».

( Читать дальше )

Если кого вдохновило сообщение smart-lab.ru/blog/680086.php, тому не обойтись без книги «NUMERICAL RECIPES. The Art of Scientific Computing. Third Edition». Качайте, пока дают

www.e-maxx-ru.1gb.ru/bookz/files/numerical_recipes.pdf
Бесплатные исходники к ней github.com/blackstonep/Numerical-Recipes
Программа svd.h из этого набора решает задачу наименьших квадратов для построения индикатора полиномиальной регрессии вместо примитивных скользящих средних.
Хорошее объяснение математической подоплёки в книге «Машинные методы математических вычислений. Форсайт, Малькольм, Моулер» en.booksee.org/book/445129
Ещё лучше — «Линейная алгебра и её применения» Гилберт Стренг
fileskachat.com/download/20151_887581203f10b39b3d7f6b84caf48a63.html
«Linear Algebra and Its Applications 4ed»
www.astronomia.edu.uy/progs/algebra/Strang- Linear_algebra_and_its_applications.pdf

Для использования программы svd.h из «NUMERICAL RECIPES» нужны тривиальные дополнения — транспонирование и перемножение матриц. Набор программ можно дополнить самодельным файлом utils.h и разместить в нём такой код:

#include <assert.h>
template <class T>
class NRdiagonal: public NRvector<T> { using NRvector<T>::NRvector; };

template <typename T>
void Multiply (const NRdiagonal<T>& a, const NRvector<T>& b
    ,NRvector<T>& c) {
  int m = a.size();
  assert (m == b.size());
  c.resize (m);
  for (int i = 0; i < m; ++i)
  c[i] = a[i] * b[i];
}
template <typename T>
void Multiply (const NRmatrix<T>& a, const NRvector<T>& b
    ,NRvector<T>& c) {
  int m = a.nrows(); int n = a.ncols();
  assert (n == b.size());
  c.resize (m);
  for (int i = 0; i < m; ++i) {
    c[i] = 0;
    for (int j = 0; j < n; ++j)
      c[i] += a[i][j] * b[j];
  }
}
template <typename T>
void Transpose (const NRmatrix<T>& a, NRmatrix<T>& b) {
  int m = a.nrows(); int n = a.ncols();
  b.resize (n, m);
  for (int i = 0; i < n; ++i)
    for (int j = 0; j < m; ++j)
      b[i][j] = a[j][i];
}
template <typename T>
void PrintVector (char* hdr, const NRvector<T>& vec) {
    cout << hdr << '\n';
  for (int i = 0; i < vec.size(); ++i)
    cout << " " << vec[i];
  cout << '\n';
}

( Читать дальше )